On the Extension of the DIRECT Algorithm to Multiple Objectives

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ On the Extension of the DIRECT Algorithm to Multiple Objectives © The Author(s) 2020 Published version Lovison, Alberto; Miettinen, Kaisa Lovison, A., & Miettinen, K. (2021). On the Extension of the DIRECT Algorithm to Multiple Objectives. Journal of Global Optimization, 79(2), 387-412. https://doi.org/10.1007/s10898- 020-00942-8 2021 Journal of Global Optimization https://doi.org/10.1007/s10898-020-00942-8 On the Extension of the DIRECT Algorithm to Multiple Objectives Alberto Lovison1 ·Kaisa Miettinen2 Received: 1 March 2019 / Accepted: 7 August 2020 © The Author(s) 2020 Abstract Deterministic global optimization algorithms like Piyavskii–Shubert, direct,ego and many more, have a recognized standing, for problems with many local optima. Although many single objective optimization algorithms have been extended to multiple objectives, completely deterministic algorithms for nonlinear problems with guarantees of convergence to global Pareto optimality are still missing. For instance, deterministic algorithms usually make use of some form of scalarization, which may lead to incomplete representations of the Pareto optimal set. Thus, all global Pareto optima may not be obtained, especially in nonconvex cases. On the other hand, algorithms attempting to produce representations of the globally Pareto optimal set are usually based on heuristics. We analyze the concept of global convergence for multiobjective optimization algorithms and propose a convergence criterion based on the Hausdorff distance in the decision space. Under this light, we consider the well-known global optimization algorithm direct, analyze the available algorithms in the literature that extend direct to multiple objectives and discuss possible alternatives. In particular, we propose a novel definition for the notion of potential Pareto optimality extending the notion of potential optimality defined in direct. We also discuss its advantages and disadvantages when compared with algorithms existing in the literature. Keywords Global convergence ·Multiobjective optimization ·Multiple criteria optimization ·DIRECT algorithm ·Deterministic optimization algorithms Mathematics Subject Classification 90C29 ·90C26 ·90C30 This paper deepens and utilizes some ideas in the extended abstract Exact extension of the DIRECT algorithm to multiple objectives [26]. This research was partly funded by the Academy of Finland (Grant no. 287496). BAlberto Lovison [email protected] Kaisa Miettinen [email protected] 1Department of Mathematics “Tullio Levi-Civita”, Università di Padova, Via Trieste 63, 35121 Padova, Italy 2University of Jyvaskyla, Faculty of Information Technology, P.O. Box 35 (Agora), FI-40014 University of Jyvaskyla, Jyvaskyla, Finland 123 Journal of Global Optimization 1 Introduction We consider the class of exact (or deterministic) global optimization algorithms and their possible extensions to the multiobjective optimization case. Among their characteristics, these algorithms exhibit appreciable convergence speed when the dimension of the decision space is moderate, but their specific strength is global convergence, the proved ability of approximating the global optimum of a black-box function with an arbitrary accuracy as e.g., in [22,30,32,35,36,40,41]. Since many real-world problems have multiple conflicting objectives to be optimized simultaneously, we concentrate here on multiobjective optimization. Such problems have so-called Pareto optima, where none of the objective functions can be improved without impairing some of the others. Some exact single objective optimization algorithms have inspired or have been used as components in multiobjective optimization algorithms like [9,16,21,24,25]. Nevertheless, the class of exact, global multiobjective optimization algorithms has not been developed to the same extent as its single objective counterpart. Indeed, in most cases, 1. asingleobjectiveexactglobalalgorithmisemployedafterthemultiobjectiveoptimization problem has been scalarized (e.g., [3,4,45]) or 2. an exact global algorithm is hybridized at some level with one or more non-deterministic principles (e.g., [4,21,42]). As a result, these algorithms can suffer from the following drawbacks: 1. If scalarization is used, a (discrete) subset of the Pareto optimal set will be accurately approximated and the quality of the approximation depends on the number of scalarized problems solved. 2. If hybridization with heuristics is used, a wide subset of the set of Pareto optima can be approximated, but heuristics cannot guarantee optimality. Moreover, in general, the convergence rate will be slower than in scalarization based algorithms. In this paper, we are interested in globally convergent deterministic multiobjective optimization algorithms. By saying that a multiobjective optimization algorithm is globally convergent, we mean that it should be able to produce an arbitrarily accurate representation of all the connected components of a globally Pareto optimal set, given a sufficiently large number of function evaluations. As mentioned above, although global convergence and deterministic characters are desirable features, it is not easy to attain them both in algorithm design. Furthermore, the question ofglobal convergence hasnot beenconsideredinmany papersin thedomain ofmultiobjective optimization [4,12,19,23–25,31,41]. The first contribution of this paper is clarifying the notion of global convergence in the case of multiple objectives on the basis of the Hausdorff distance between sets. Indeed, the use of set-based distances like Hausdorff is not relevant for single objective optimization algorithms, because the global optimum for single objective optimization is a unique value in the objective space, corresponding to a set of minimizers in the decision space that in a typical case consists of a unique point, but can consist of more than one point or even an infinite number of points in pathological cases. Therefore, the Euclidean distance between the exact minimizer and an approximated minimizer is the only necessary concept for establishing a convergence speed. This is not the case in multiobjective optimization, where in the generic case, the set of Pareto optima is a (k−1)-dimensional submanifold of the objective space, if k>1isthe 123 Journal of Global Optimization number of objectives. To establish the convergence speed of a multiobjective optimization algorithm, it is necessary to measure the distance between the exact Pareto optimal set and the approximated Pareto optimal set. This issue is not addressed univocally in the literature. The second contribution is an analysis and a discussion of the well-known direct algorithm [20] and its extensions to multiple objectives. We focus on direct because it has been recognized as an efficient algorithm for medium sized problems in the decision space. Indeed, direct doesnotsuffer muchofthe curseof dimensionalityunlike, e.g.,thePiyavskii–Shubert algorithm[33,37].Furthermore,direct doesnotrequireanyglobalinformation(e.g.,aglobal Lipschitz constant) to be globally convergent and, finally, it has been a starting point for many efficient single objective optimization algorithms. In the original direct algorithm, the feasible set is the d-dimensional hyperinterval D= [0,1]d⊆Rd.Dis covered by small hyperintervals D1,...,DN,whereD=[a() 1,b() 1]× ···×[a() d,b() d],with 0 ≤a() i<b() i≤1,1≤i≤d, such that ∪N =1D=D.direct contains two fundamental steps: a subdivision step and a selection step. In the division step, one of the small hyperintervals Dof the covering of Dis partitioned into smaller hyperintervals. In the selection step, the hyperintervals that are likely to contain the minimum are selected for further subdivision. The selection rule is based on a notion of potential optimality. This notion cannot be extended to multiple objectives in a straightforward way. We present, therefore, a review and a discussion of the existing extensions of direct.Next, we propose a novel formulation for the notion of potential Pareto optimality, which is as faithful as possible to the original single objective version. We also discuss the strengths and weaknesses in comparison with the existing alternatives. As far as the rest of this paper is concerned, we provide background in Sect. 2by introducing the main concepts and notations used in the rest of the paper. Furthermore, we give a literature survey of deterministic, global multiobjective optimization algorithms, in particular, those based on direct. For the convenience of the reader, we also briefly resume the single objective direct algorithm. In Sect. 3, we discuss the notion of global convergence in single objective optimization, as it has been interpreted in the literature. Then we extend the focus to how it should be formulated to guarantee that an algorithm will accurately represent every possible connected component of a Pareto optimal set in multiobjective optimization. In Sect. 4, we discuss how the principles of direct have been extended to multiple objectives in the literature and analyze possible alternatives. These ideas are collected and presented in the proposed multidirect algorithm. In Sect. 6, we discuss the behavior of multidirect,the multiobjective extension of the Piyavskii–Shubert algorithm mPS [25] and another direct extension developed recently [1] on a simple but paradigmatic example defined in [19]. Finally, we conclude in Sect. 7. 2 Background 2.1 Some concepts of single and multiobjective optimization There is a fundamental aspect that distinguishes the global optimization of a single function and several objective functions. Let us first consider the optimization problem with a single objective continuous function f:D→R,D=[0,1]d⊆Rd,minimizex∈Df(x), (1) 123 Journal of Global Optimization where Dis called a feasible set in the decision space Rdand there are no other constraints. Because of compactness, the set of globally optimal objective function values for problem (1) is a single value that often corresponds to a unique optimal point in the decision space. In some cases, even an infinite number of optimal points can correspond to the unique optimal value. However, this situation can be proved to be nongeneric, i.e., for a dense set in a suitable set of functions, there exists only one optimal value with a unique corresponding optimal point in the decision space. This is a consequence of Sard’s theorem and Morse’s lemma [17]). Therefore, the global convergence of single objective optimization algorithms is established by requiring that the globally optimal value is approximated as accurately as desired, given an infinite number of function evaluations. In the generic case, this determines a unique optimum also in the feasible set D. In the case of multiple continuous objectives, i.e., f:D→Rk,D=[0,1]d⊆Rd,minimize x∈Df(x), (MOP) where f(x)=(f1(x),... fk(x)) is called an objective vector in the objective space Rk,in general, it is not possible to completely order multidimensional vectors. Indeed, one vector can be better than another one when considering their first components, but the converse may happen when comparing the second components. Therefore, none of the two vectors can be said to be better than the other one. Only in special cases, say, pathological ones, there could exist a point which is the global optimum with respect to all the objective functions at the same time. In such cases, multiobjective optimization methods are not needed. To dealmorerigorouslywithoptimality,itisnecessaryto considerthenotionofdominance [18,29]. An objective vector zin Rkdominates another objective vector w∈Rk,andwe write z≺w,ifzi≤wifor all iand there exists jsuch that zj<w j. We say that a point xisa(global) Pareto optimum if f(x)is nondominated by any other f(y),foryin the feasible set D. The set of Pareto optima form a Pareto set P.Thesetof corresponding objective vectors, i.e., {f(x):xis a Pareto optimum}is called a Pareto front. If a point xis nondominated by any other point in one of its neighborhoods, it is called a local Pareto optimum. If a Pareto optimum is not a minimizer for one of the objective functions considered individually, we call it a tradeoff point. If objective functions have been evaluated only for a finite number of points S={x1,...,xN}⊆D, we can consider the subset of Sconsisting of the Pareto optimal points only with respect to Sand not with respect to the entire D. These points cannot be guaranteed to be Pareto optimal, although they are used as an approximation of the Pareto set. We call this subset a nondominated subset of Sand denote it by ND(S). We recall also a widely used concept related to Pareto optimality [29, Definition 2.5.1]. A point x∈Dand the corresponding objective vector f(x)are weakly Pareto optimal if there does not exist another point y∈Dsuch that fj(y)< fj(x)for all j=1,...,k. 2.2 Literature survey In what follows, we discuss multiobjective optimization algorithms currently available in the literature that are either deterministic or related to the problem of finding the set of global Pareto optima. Special attention is dedicated to algorithms extending or involving direct. NormalBoundary Intersection(nbi)[10]isa deterministicalgorithm capableof producing an evenly spaced discretization of the Pareto front by defining a family of scalar subproblems of the original multiobjective optimization problem. nbi is based on continuation, i.e., it finds connections among different Pareto optima. nbi starts from a Pareto optimal point, which is 123 Journal of Global Optimization (b)(a) Fig. 1 Panel a: application of sicon [23]toL&H2x2 [19]. The algorithm detects correctly all the different connected components of the Pareto set. Panel b: application of the NBI algorithm to the L&H2x2 problem. The starting set is a regular grid. The algorithm correctly traces only one of the connected components of the Pareto optimal set, missing the second component. (Numerical experiment on modeFRONTIER® courtesy by E. Rigoni at ESTECO s.r.l.). the solution of one of the scalar subproblems, and identifies the local connected component of the Pareto set containing the starting point. In recent years NBI has been extended and generalized in other algorithms, see for instance [3]. As a drawback, in nonconvex cases, nbi possibly misses multiple isolated components, giving rise to an incomplete representation of the Pareto set. This issue is observed in [10] and discussed in [23], where a continuation algorithm detecting correctly all the separated connected components, called sicon, is introduced. See Fig. 1a for an illustration of such a phenomenon for a nonlinear and nonconvex problem, called L&H2x2 [19] (in both decision and objective spaces). In it, the Pareto set consists of two components, although their images are mapped in a single connected Pareto front. Another deterministic continuation algorithm is proposed in [34] solving the first order Karush–Kuhn–Tucker (KKT) conditions for Pareto optimality. Points satisfying the KKT conditions are called substationary. Substationary points found are approximated by a collection of boxes obtained by an iterative subdivision of the feasible set of the problem. This algorithm has two main drawbacks: first it approximates only one connected component of the Pareto set and, second, the set of substationary points is a strict superset of the Pareto set. In [12], the authors extend the ideas of [34] and propose an algorithm that, by using the properties of a suitably defined dynamical system, approximates a set of substationary points (as defined above). The output of the algorithm is a covering of a set composed by a collection of boxes. The authors discuss convergence towards a global Pareto set in a similar way as done in this paper by adopting the notion of a Hausdorff distance between sets. The main difference to our algorithm is that for the global convergence, the connectedness of the Pareto set is required and, therefore, multiple connected components of the Pareto set cannot be handled properly. Singular Continuation (sicon) proposed in [23] is another continuation algorithm based on both the first and second order optimality conditions proposed by Smale in [38]. It has some similarities with [12]. To be more specific, sicon aims at approximating the set of substationary points. In addition, sicon takes into consideration the second derivatives and, therefore, substationary points can be separated in the multiobjective analogue of minima, maxima and saddle points. When compared with the previous strategies, sicon can produce a piecewise linear approximation [2] of the Pareto set in the form of a triangulated surface, 123 Journal of Global Optimization with a provably Newton-like superlinear accuracy. Furthermore, the algorithm allows for determining all the connected components of the Pareto set, provided that the starting sampling is dense enough. Indeed, in Fig. 1b, two connected components are detected, but if the sample of points had been sparser, the algorithm could have missed the smaller component. The sicon algorithm relies on a Delaunay tessellation of the feasible set and can suffer from the curse of dimensionality. Because of how it is defined, it is not possible to add new points and refine the structure automatically, although extensions in this direction are possible. In [24,25], the question of global convergence is discussed on the basis of the Hausdorff distance in the decision space. In [24], an algorithm globally converging towards the substationary set is proposed, while in [25] an extension of the Piyavskii–Shubert globally convergent algorithm is introduced and the convergence towards the global Pareto set is proven assuming that a global Lipschitz constant is known for each objective function. Lipschitz optimization for univariate problems is considered for the bi-objective case in [41, Chapter 7], in [46,47], and for the multivariate case in [31]. In [19], the authors discuss the problem of correctly identifying different components of the Pareto set that may superimpose when nonconvex objective functions are involved. The aim is to produce globally correct representations of the Pareto set by means of triangulated surfaces, extending a previous nonglobal algorithm [18]. A deterministic approach making use of the Hausdorff distance in the objective space is called the Non Uniform Space Covering Method [14,15]. It provides an ε-approximation of the Pareto front in the sense that the Hausdorff distance between the Pareto front and its approximation is guaranteed to be smaller than ε>0, given some suitable assumptions on the estimation of lower bounds. MultiGLODS [9]isa hybrid local-global deterministic algorithm combining Direct Multisearch [7], a derivativefree multiobjective local optimization algorithm and GLODS [8], a space-filling optimization method that manages efficiently multiple sequences possibly converging towards the same Pareto optimal point. Nevertheless, the authors focused on local convergence features and proved local convergence towards at least one Pareto optimal point. Let us finally survey some algorithms related to the direct algorithm, although some of them make partial use of heuristics and are not, thus, fully in our scope of studying deterministicalgorithms. In[1,45],theauthorspropose threealgorithms calledMO-DIRECT, derived from direct by employing either the ranks of the nondominated sorting [11], the hypervolume indicator or the nondominance (ND) concept in the objective space augmented by an extra dimension consisting of the hyperinterval radius. Several numerical performance tests have been implemented and used for measuring the capabilities of these algorithms. For the purposes of comparison with the algorithm proposed here, we notice that only ND does not make use of scalarizations. A proof of global convergence in the sense that the algorithm will produce an everywhere dense sampling of the feasible set is provided. The code is available and we provide a theoretical and numerical comparison with MO-DIRECT-ND in the subsequent sections. A series of papers [5,6] discusses applications of an extension of direct to damage identification problems with a numerical comparison with MO-DIRECT [1,45]. The algorithm proposed uses a scalarization of the multiobjective optimization problem by assigning to a vector vin the objective space the rank R(v) given by the nondominated sorting. The papers are application-oriented and the problem of global convergence is not considered. A multiobjective extension of direct (called modir) is proposed in [4] exhibiting an engineeringapplication.The algorithmuses a non-explicit scalarizationofobjectives because it selects as potentially optimal those hyperintervals that are potentially optimal according to at least one of the objectives. This definition excludes the selection of tradeoff points. Nevertheless, the authors provide a proof of global convergence in the sense of everywhere 123 Journal of Global Optimization dense sampling in an infinite time. As an alternative formulation, the authors also propose a selection rule based on the nondominated set in the augmented feasible set using the hyperinterval radius as in MO-DIRECT-ND. The authors also propose a hybridization with a stochastic local search strategy called DFMA. An evolutionary hybrid algorithm called NSDIRECT-GA is proposed in [42]. There, direct is combined with MOGA and compared with NSGA (both evolutionary algorithms, see [11]). 2.3 The DIRECT algorithm direct [20] is quite a natural generalization of the Piyavskii–Shubert algorithm [33,37]and outperforms its ancestor in several senses. First, direct does not require the existence of a global Lipschitz constant, which usually is not known, and second, direct exhibits a faster convergence. In the following subsections, we outline the working principles of direct. This will make more comprehensible the ideas behind the multidirect algorithm to be proposed. First, we consider a problem with a Lipschtiz continuous single objective function ffor which the Lipschitz constant is unknown. 2.3.1 Selection rule In direct, the feasible set Dis the d-dimensional hyperinterval D=[0,1]d⊆Rd, covered by a collection of hyperintervals D1,D2,...,whereDj=[a(j) 1,b(j) 1]×···×[a(j) d,b(j) d], with 0 ≤a(j) i<b(j) i≤1,1≤i≤d. The hyperintervals D1,D2... can overlap only along the faces. The radius σjof the jth hyperinterval Djis defined as σj=maxi=1,...,d b(j) i−a(j) i 2. Some of these hyperintervals are selected for subdividivision on the basis of a Lipschitz-like criterion. Assuming that the function fhas been evaluated at the point in the center of the jth hyperinterval, it is possible to estimate a lower bound for fby fixing a rate of change K >0 (which has a role similar to a Lipschitz constant) by setting j=f(xj)−Kσj,wherexjis the point at the center of the jth hyperinterval and σjis the radius of the hyperinterval. By ordering the sequence of jwe obtain a ranking of the hyperintervals that are most likely to contain the global optimum of the function f. This ranking can vary when different rates of change Kare chosen. If Kis small, hyperintervals with low values of f(xj)will be highly ranked, even if they have small σj. On the other hand, when Kis large, hyperintervals with a large radius σjwill be highly ranked. We call potentially optimal a hyperinterval jwhose lower bound is optimal for a certain fixed rate of change K, i.e., f(xj)−Kσj≤f(xi)−Kσifor all i.(2) We can show that all potentially optimal hyperintervals correspond to points in the lower face of the convex hull of the following set of points as the rate of change Kvaries from 0 to ∞, as depicted in Fig. 2: σj,f(xj)j=1,2,...⊆R2.(3) We call the plane {(σ, f):σ, f∈R}as an augmented plane of values. At each iteration, direct picks these potentially optimal hyperintervals on the convex hull and promotes them for further subdivision in smaller hyperintervals. It is typical that several hyperintervals are 123 Journal of Global Optimization (a) (b) Fig. 2 direct selection rule. aFunction values evaluated at the center of a hyperinterval versus hyperinterval radius σ=maxi=1,...,d(bi−ai)/2. The lower bound corresponding to a fixed rate of change Kis the intersection of the line y=f(cj)+K(σ −σj)with the yaxis. bHyperintervals corresponding to points lying on the lower, right-hand side of the convex hull are potentially optimal to be subdivided, because as Kvaries, the lower bounds of each hyperinterval vary and, therefore, the ranking changes. This selection rule indeed selects a subset of the Pareto set of a biobjective optimization problem, where the objective functions are the hyperinterval radius, to be maximized, and the function value, to be minimized. The subset in fact corresponds to the convex hull of the Pareto front for this problem. For practical reasons, a further rule of selection is adopted: the lower bound computed should improve the current lowest value computed by at least a fixed quantity , i.e., f(xj)−Kσj≤fmin −|fmin|.(4) This extra rule prevents oversubdivision of the smallest hyperintervals which could easily lead to the depletion of the available computational resources without significantly improving the algorithm performance. 2.3.2 Division rule A crucial aspect of the algorithm is the way hyperintervals are divided in smaller ones. The original algorithm proposes first to compute two extra test points along each coordinate direction of the feasible set, and then on the basis of the values computed, the hyperinterval is divided in thirds of equal size by first considering the directions where function values are more promising. So, if we call xi1and xi2the new points considered along the ith direction and fi1and fi2the corresponding function values, directions are ranked from the lowest to the highest min {fi1,fi2}. Then, the most promising directions according to this criterion are divided in thirds at first, such that the best values are contained in the subintervals with the largest diameter. This feature improves the convergence speed, because best rated directions are more likely to be subdivided earlier which increases the density in the neighborhoods of test points with highest ranking. The division rule is exemplified in Fig. 3. A further detail is that actually not all the directions are divided but only the ones corresponding to sides with the largest length. We notice that because of the division in thirds at 123 Journal of Global Optimization Fig. 6 Potential Pareto optimality for three hyperintervals on different σ-levels. aVector lower bounds for H1 and H3both dominate the vector lower bound for H2for rates of changes α1and α2with the same magnitude. bIf α1α2, the vector lower bound for H2is no more dominated. cCorrespondingly, in the augmented space we see that some of the dominance rays intersect the higher σ-levels in points that are nondominated in the σ-level. dIf there are two objective vectors F3and F4in the upper level dominating the intersection of the dominance cone of H2with their σ-level, then H2is not potentially Pareto optimal. Notice that the most important rays (in blue) are at the extrema of the virtual Pareto front (i.e., the red line). Colors are available in the online version of the paper σ→  F(A,σ),σ= Fi(σ), σ :=  F+Fi− F σi−ˆσ(σ −ˆσ),σ.(7) Note that  Fi(ˆσ) = Fand  Fi(σi)=Fi, as claimed. Moreover  Fi(0)is the vector lower bound both for  Fand Fiwhen Aas above is used. Proposition 1 A hyperinterval  H is potentially Pareto optimal if and only if there exists a set of rates of change A =(α1,...,α k)such that the intersection of the ray σ→ (FA(σ), σ ) := ( F+A(σ −ˆσ),σ) with every σ-level is nondominated by the objective vectors Ficorresponding to that σ-level, i.e., such that σi=σ. Proof We observe that  His potentially Pareto optimal if and only if there exists A such that  FAis nondominated in the set of lower bounds associated with A, i.e., FA i=Fi−Aσi|i=1,2,.... Moreover,  FA≺FA i⇐⇒  F−Aˆσ≺Fi−Aσi⇐⇒  F+A(σi−ˆσ) ≺Fi.(8) 123 Journal of Global Optimization But  F+A(σi−ˆσ)is exactly the intersection with the σi-level of the line σ→  F(A,σ),σ passing through ( F,ˆσ). The claim follows immediately.  We introduce a subset of the objective space very important for the geometric characterization of potential Pareto optimality. Definition 7 Let E=v∈Rk Fj(0)v Ffor some j∈{1,...,t}. We refer to the following set as a virtual Pareto front: VP[ F]:=(v, 0)∈Rk+1v∈E,vis weakly Pareto optimal in E.(9) The cone of lines connecting ( F,σ) with the virtual Pareto front is called dominance cone of  H. In Figs. 5c, 6c, d, and 7b, the virtual Pareto front is highlighted in red, while the dominance cone is the green surface (colors available in the online version of the paper). Proposition 1requires that a high-dimensional infinite set of possible rates of change A=(α1,...,α k),αi≥0, is checked in order to determine whether  His potentially Pareto optimal. Nevertheless, it suffices to check only for Ain a lower-dimensional subset. Proposition 2 A hyperinterval  H is potentially Pareto optimal if and only if there exists A=(α1,...,α k)such that (A,ˆσ)is in the dominance cone of  H, and the intersection of the ray  FA(σ) :=  F+A(σ −ˆσ) with every σ-level is nondominated in the σ-level, for every level σ>ˆσ. Proof If such a ray A∈Rkexisted, it would be possible to find one  Ain the dominance cone by picking the ray passing by the intersection of the virtual Pareto front and the line connecting  FAand  F. Indeed, consider λAin place of A, with 0 <λ<1. By continuity, there exists a neighborhood of 1 such that  FλAis still nondominated, as well as for the intersections of the dominance ray with the σ-levels with σ>ˆσ. The infimum of such values of λ, by definition, corresponds to a point in the virtual Pareto front.  Proposition 1is the multiobjective extension of the geometric criterion in (2). As equation (2) requires an infinite number of checks, i.e., all values α,0<α<∞, also Propositions 1 and 2would require to test an infinite number of A’s, at least to say that  His not potentially Pareto optimal. In the single objective case, it is possible to limit oneself to the slopes of the lines connecting points in the convex hull of the augmented function values, therefore to reduce to a finite number of tests. This is possible also in the multiobjective case. Definition 8 Tipping vectors T ip are elements in the virtual Pareto front which are nondominated if the dominance relation is reversed by using in place of ≺, i.e., Tip = −ND(−VP[ F]). Remark 2 In the biobjective case, tipping vectors are the concave vertices of the virtual Pareto front, plus the two extrema. In Figs. 5,6,7, the tipping vectors and the corresponding dominance rays are highlighted in blue (colors are available in the online version of this paper). Solid blue lines correspond to intersections that are nondominated in all the higher σ-levels, while dashed lines correspond to dominated intersections in some of the higher σ-levels. Theorem 1 An hyperinterval  H is potentially Pareto optimal if and only if there exists a tipping vector T such that the ray connecting ( F,σ) with T is not intersecting the regions dominated by the remaining vectors Fifor each of the σ-levels with σ>σ. 123 Journal of Global Optimization Fig. 7 Geometric characterization of potential Pareto optimality for six hyperintervals. Here, F1and F2belong to the same σ-level, as well as F4,F5and F6, i.e., σ1=σ2<σ 3<σ 4=σ5=σ6. Since F1and F2are nondominated, they are potentially Pareto optimal. While F4,F5and F6are nondominated in the higher σ-level, they are also potentially Pareto optimal, according to Remark 1. The problematic hyperinterval is H3. By inspecting the tipping vectors (blue dots) and the respective dominance rays we see that the central ray intersects the higher σ-level in a nondominated point. Hence, also H3is potentially Pareto optimal Proof For every vector  F−Aˆσin the virtual Pareto front, there is (at least) a tipping vector Tsuch that  F−Aˆσ≺T. Concerning the intersections of the corresponding dominance cones with higher σ-levels, they are FA σ= F+A(σ −ˆσ) and Tσ= F+(T− F)σ−ˆσ 0−ˆσ. Explicit calculation with σ>ˆσand  F−Aˆσ≺Tshows that Tσ≺FA σand, therefore, if FA σis nondominated in the σ-level, also Tσis nondominated.  Example 5 In Example 4, all hyperintervals are potentially Pareto optimal. This appears also in Fig. 6c, where we observe that the dominance rays corresponding to two tipping vectors are nondominated in the higher σ-level (solid blue lines; colors available online). Consider four hyperintervals H1,H2,H3,H4, with objective values F1,F2,F3,F4∈R2, with F1≺F2, F2≺F3and F2≺F4as in Fig. 6d. In this case, H2is not potentially Pareto optimal. Indeed, all the dominance rays passing by vectors in the virtual Pareto front are dominated in the higher σ-level (dashed blue lines; colors available online). Example 6 Consider six hyperintervals H1,H2,H3,H4,H5,H6, with F1,F2≺F3,F3≺F5 as in Fig. 7, with σ1=σ2<σ 3<σ 4=σ5=σ6. In this case, H3is potentially Pareto optimal, because the dominance ray (solid blue line) corresponding to the central tipping vector does not intersect the dominated region in the higher σ-level. The dominance rays of the remaining tipping vectors are dominated in the higher σ-level (dashed blue lines). 4.3 Supplementary condition for potential Pareto optimality An important condition in the definition of potential Pareto optimality is (4), which depends on a small parameter . However, according to the authors of direct ([20]), the algorithm seemsquiteinsensitiveonthe magnitudeoftheparameter,provideditissmallwhencompared to the magnitudes of the ranges of the objective function values. This parameter dependent condition is critical for avoiding that the algorithm wastes excessive computational resources 123 Journal of Global Optimization on exploring regions with very small potential improvement. This condition, therefore, is important for efficiency. As far as we know, in the literature, condition (4)ofdirect has never been adapted and implemented for multiobjective cases. If we denote by |F|the vector of absolute values of the components of F, i.e., |F|:= (|f1(c)|,...,|fk(c)|), we require the vector lower bound  FAbeing nondominated by the objective vectors Fidiminished by a quantity proportional to and |Fi|:  FAis nondominated by the set Fi−|Fi|i=1,...,N := (f1(ci)−|f1(ci)|,..., fk(ci)−|fk(ci)|,)i=1,...,N.(10) Operatively, this has no effect on the vectors belonging to the σ-level with maximal σ. For the vectors belonging to the intermediate σ-levels, it is necessary to substitute the tipping vectors Tin Theorem 1with T−→ T+|T|.(11) This does not implement exactly condition (10), but it appears as a valid approximation. Notably, condition (10) implies that also nondominated vectors belonging to the lowest σlevel, which were automatically selected as potentially Pareto optimal, should also be filtered. For these points, however, there is no tipping vector defined, because there are no dominating vectors from which to compute a virtual Pareto front. Therefore, we use the objective vector  Fas a tipping vector, then we proceed as above by subtracting a vector proportional to |F| and , as for the previous case. The dominance ray is the line passing through ( F,σ) and ( F−|F|,0). If the dominance ray intersects all the higher σ-levels in a vector that is nondominated for all possible σ, then the hyperinterval  His potentially Pareto optimal. As we will see in Sect. 6, this condition becomes more and more influential in the later iterations. There, a large part of the points on the Pareto front, which becomes denser and denser, is excluded from selection because the points do not bring a relevant potential improvement. 4.4 Discussion and comparison with the selection rules of other algorithms Alternativeextensions ofthe concept ofpotential Paretooptimality are consideredin [1,4,26]. In [4], a hyperinterval His deemed potentially Pareto optimal if it is potentially optimal in the scalar sense (2)forat least one of the objective functions. In [26], an exact extension of the geometrical interpretation of potential optimality is consideredasaconvexhullofthesetofaugmentedobjectivevectors.In[1],thenondominated set of the set of the augmented objective vectors is used as the set of potential Pareto optima. The same definition is considered also in [4] as a possible alternative. By considering Example 6, we notice that for the algorithm of [4], the potential Pareto optima are only H1,H2,H4and H6. Here, although H5is a nondominated point in the σlevel, it is not deemed potentially Pareto optimal because it is not optimal for at least one of the objectives. In the generic case, the subset of the Pareto set composed by tradeoff points is a (k−1)-dimensional submanifold of the decision space, while points that are optimal for at least one objective considered individually, form a discrete and finite set. Therefore, the definition in [4] excludes a large part of potentially relevant hyperintervals. Also according to [26], H1,H2,H4and H6are deemed potentially Pareto optimal. In this case, H5is excluded because it belongs to a nonconvex part of the nondominated set of the σ-level. Nevertheless, the notion in [26] is less restrictive because it does not exclude at least 123 Journal of Global Optimization (a) (b) Fig. 8 In direct, a ranking among directions is provided by objective function values (a). In the multiobjective case, a total ordering is not available a priori (b) the convex parts of the nondominated sets, while [4] can pick as potentially Pareto optimal only the extrema of the subfronts. On the other hand, if we consider the nondominated subset of the augmented objective vectorsasin[1], we can say that all hyperintervals H1,...,H6are potentially Pareto optimal, which is in agreement with Definition 5. However, Definition 5excludes H2in Example 4, while for [1]allH1,...,H4are potentially Pareto optimal. We can conclude that, when comparing the new definition 5with the definitions used in the previous algorithms, there is a possibility that [4]and[26] result in slower in convergence, because they exclude promising hyperintervals in early iterations. The algorithm of [1] can be slower in convergence, because it could waste some computational resources for hyperintervals that are less likely to contain Paretooptimalpointsthanotherhyperintervals.Theeffectivenessofthesedifferentdefinitions of potential Pareto optimality are compared in Sect. 6. 4.5 Division rule The challenges of extending the division rule to multiple objectives are related to the fact that in the original division rule for the single objective case, a total ordering criterion among the directions in the decision space is necessary. This total ordering is obtained by testing values computed at both sides of the central point, as illustrated in Figs. 3and 8a. Structurally, this is a feature impossible to be reproduced in multiobjective optimization, because test values computed along different directions are vector-valued. Thus, they cannot be compared univocally. This is exhibited in Fig. 8b, where it is clear that comparisons are not always possible. The central point, for instance, can be compared with points in the first and third quadrant of the plane, but not with points in the second and fourth quadrant. In the literature, some criteria have been proposed to extend the division rule using, for instance, the Euclidean distance of the test point values from the objective functions value of the center point of the ancestor hyperinterval [1,45]. Other choices are possible, e.g., 123 Journal of Global Optimization using an achievement scalarizing function [43] with the central point value as a reference point. Alternatively, the hypervolume [13] corresponding to each new test point can be used, producing a scalar value useful for ranking. These algorithms, however, make use of a more or less explicit reference to the individual scales at which the objective functions are represented. The Euclidean distance and the hypervolume depend on the function scaling and the achievement scalarizing function explicitly requires that some weights are provided to combine the functions. Instead, a lexicographic ordering could be used, but also in this case an explicit, and arbitrary, ranking among objectives is introduced, which shares the same drawback observed above. Of course, in a general case, an a priori ranking among the objective functions, function rangesor respectiveglobal Lipschitzconstants isnot provided anddoesnotexistinauniversal sense. This means that each of the above-mentioned rankings is somehow arbitrary, i.e., they introduce a bias that in some cases may be negligible, but in principle it can impair the functioning and the convergence features of the algorithm. A possible neutral approach could be an iterative estimation of the Lipschitz constants (or a reciprocal scaling) on the basis of the values computed so far, and an attractive choice could be that the directions showing the largest estimated Lipschitz constant could be divided earlierthantheremainingdirections tofollowpossibleemerginganisotropyinthefeasible set. Nevertheless, it is preferable, and conceptually cleaner, that the sequence of newly computed points is completely independent of any ranking or scaling of objective functions. This is in agreement with the issues discussed in the introduction regarding the existing algorithms extending or making use of direct and other exact global optimization algorithms. A possible approach avoiding a total ordering between objectives is to use nondominated sorting [39], which is a partial ordering. In nondominated sorting, we extract from a finite discrete set of vector test values the nondominated ones, to form a so-called first subfront. After removing the first subfront from the starting set we select again the nondominated vectors from the remaining vectors to form a second subfront, and so on until the original set is completely emptied. Test values belonging to the same subfront are given the same rank. The computation of each subfront relies only on Pareto dominance and, therefore, it does not depend on the function scalings or on any other arbitrary ranking between the objective functions. Nondominated sorting is illustrated in Fig. 9a, where the first subfront is the blue one, the second is the orange, etc. By adopting such a partial ordering among the directions in the feasible space to be subdivided, we modify the original division rule by splitting at the same time directions belonging to the same subfront, i.e., having the same rank. This requires that some extra test point is computed at the center of the extra hyperintervals generated in this way, as can be seen in Fig. 9c. In addition to avoiding the use of any ranking between directions or scaling between objectives, this partial ordering avoids introducing complicated or arbitrary adaptive rules. On the other hand, it involves some computational burden, that has to be taken into account. The computational cost of extra points may be acceptable because of the following reasons. First of all, in typical applications, the dimension of the feasible space, which corresponds to the number of directions to be divided, cannot be very large, because an exhaustive search requiring global convergence, i.e., an everywhere dense sampling in an infinite time, cannot be efficiently performed in high-dimensional decision spaces. Rarely we are considering decision spaces with dimensions larger than five. Secondly, not all directions in the decision space are candidates for division but only the subset corresponding to the longest edges of the hyperinterval, as in the original direct 123 Journal of Global Optimization (a) (b) (c) Fig. 9 Division rule. aTest values are grouped in subfronts according to nondominance. bSometimes the partial ordering is sufficient for establishing a total ranking among input directions and the classical division rule can be applied. Numbers represent objective function values. Test points corresponding to the same direction belong to the same subfront and obtain the same rank. cWhen a total ranking is not available, a complete division rule is adopted. Test points corresponding to different directions belong to the same subfront division rule. Finally, nondominated sorting produces a total ordering among directions in a considerable number of cases, as can be seen in Fig. 9b. 5multiDIRECT algorithm Now we can collect and summarize the ideas discussed in the previous sections about the possible issues of extending the direct algorithm for multiobjective optimization problems. Based on the discussion we can design an algorithm to be called multidirect,whichis presented in Algorithm 1. Algorithm 1 multidirect Input: Let fbe a vector-valued function as in (MOP). Let nIter be the number of iterations to be performed. Let >0 be a small parameter. 1: Initialize the set of hyperintervals covering the domain D=[0,1]dby S:= {H}with H:= (c,F,σ) = ((0.5,...,0.5), F(0.5,...,0.5), (1,...,1)). 2: for iter =1,iter ≤ nIter ,iter ++,do 3: Group the hyperintervals according to their radius in σ-levels σL:= H=(c,F,s)s=σ. 4: For each σ-level, extract the nondominated set ND(σ L). 5: Initialize the set of candidate hyperintervals Cand with the nondominated set of the σ-level with σ= σmax, i.e., Cand := ND(σmax L). 6: In the nondominated set of σ-levels with a nonmaximal σ, add the potentially Pareto optimal hyperintervals according to Theorem 1and condition (11)toCand. 7: For every hyperinterval in Cand, select the input directions with a maximal radius, divide in thirds along each of these directions and define new test points at the center of each third. 8: Evaluate objective functions on the test points. 9: Determine the partial ordering among directions according to the nondominated sorting. 10: Apply the division rule described in Section 4.5 evaluating new test points where necessary. 11: Updatethe list Sby removing the candidate hyperintervals and addingthe newlycomputedsubintervals. 12: end for 13: Extract the nondominated set of the hyperintervals Sas the approximation of the Pareto optimal set. return ND(S). 123 Journal of Global Optimization We notice that the algorithm is globally convergent. Theorem 2 The sequence of points generated by multidirect for any function f :D→Rk becomes infinitely dense everywhere in D in an infinite time. Proof By recalling Remark 1, we note that for every iteration, the nondominated vectors in the σ-level with a maximal σare regarded potentially Pareto optimal. Therefore, in a finite number of iterations, the maximal σ-level will be emptied, i.e., σmax →0 in a infinite time. Then, also the radius of a ball contained in the feasible set and not containing any of the sampled points, i.e., the centers of the hyperintervals, is tending to zero, because this radius is smaller than the diameter of the largest hyperinterval multiplied by √d. Therefore, the sampling consisting in the centers of the hyperintervals is deemed to become everywhere dense.  6 Some numerical comparisons In this section, we demonstrate the numerical behaviour of the proposed multidirect algorithm and compare its performance to other deterministic algorithms. In the sequence of visualizations in Fig. 10, we present the application of the multidirect algorithm on a biobjective problem with two variables named L&H2×2. The problem is defined in [25] and its Pareto set and front are illustrated in Fig. 1a. We have used the output of the sicon algorithm as a surrogate for the exact representation of thePareto set, because sicon produces a set-wise approximation with a quadratic precision. Furthermore, it was possible to produce a grid of points of desired density in the Pareto set. On the other hand, we did not compare sicon with the remaining methods because it is not of the same class as multidirect. Indeed sicon makes use of first and second derivatives and it is not designed for being globally convergent. We have decided to compare multidirect with two other algorithms because they were the only globally convergent and completely fully deterministic ones available. In addition, we have also tested a version of multidirect with a division rule borrowed from one of the other algorithms. More precisely, we consider: 1. MO-DIRECT-ND with total ordering in the division rule as in the original paper [1] (green line), 2. multiPS with the Lipschitz constants computed from the analytical definition of the functions (orange line), 3. multidirect with thesame totalorderingforthedivision ruleinMO-DIRECT-ND (yellow line), 4. multidirect with partial ordering as in Algorithm 1and =0 (blue line). To compare the performance of the four algorithms, we have measured at each iteration the Hausdorff distance between the Pareto set approximation produced by the algorithm under examination and the surrogate exact Pareto set produced by sicon.InFig.11,weplotthe Hausdorff distances versus the number of function evaluations used to produce the current approximation. In Fig. 11, we also report the mesh size used in sicon, the red horizontal line, to represent the best approximation level reachable using this surrogate for the Pareto set. As can be observed, the algorithms using the partial ordering require initially a larger number of function evaluations but demonstrate better performance in the last iterations. By observing the two versions of multidirect, i.e., the yellow and blue lines in the online paper, there are no differences in performance for the first iterations. Actually the blue line 123 Journal of Global Optimization Fig. 10 Application of the multidirect algorithm to L&H2x2[19]. For every subpicture, the left panel depicts the decision space while the right panel represents the objective space. The last subpicture is a representation of the local Pareto set obtained by sicon [23] is hidden behind the yellow line. When the selection rule based on the potential Pareto optimality starts to filter some of the hyperintervals accepted by the ND rule, the number of function evaluations for every iteration gets smaller while the performances do not deteriorate sensibly. This results in with an overall better efficiency. Although the partial ordering division rule used by multidirect appears more costly in the first iterations, it leads to smaller distances from the global Pareto set in later iterations. ThemultiPS algorithmappears overall moreefficientthantheotheralgorithms,at leastwith 103to 3×104function evaluations, but it must be taken into account that global information of the analytical Lipschitz constant has been used and it is usually not available. Using larger constants implies wasting computational resources in unpromising regions, while smaller Lipschitz constants may miss to detect some non-negligible portion of the Pareto set. 123 Journal of Global Optimization Fig. 11 Global convergence rates for the extensions of the direct algorithm. Green line: MO-DIRECT- ND. Orange line: multiPS algorithm with analytical estimate of the global Lipschitz constant. Yellow line: multidirect with total ordering as in MO-DIRECT-ND. Blue line: multidirect as in Algorithm 1and =0. The red horizontal line represents the mesh size with which the exact Pareto set has been approximated, therefore the Hausdorff distance cannot go under this value We can eventually conclude that the overall best performance has been obtained by the multidirect method, although it required a considerable number of function evaluations. In particular, it seems that the selection rule based on the new notion of potential Pareto optimality is a crucial aspect in the efficiency of the algorithm proposed. 7 Conclusions We have provided a conceptual and theoretical investigation on the notion of global convergence in multiobjective optimization, with a particular emphasis on the notion of potential Pareto optimality. It is a crucial notion involved in extending the global single objective optimization algorithm direct to multiple objectives. We have discussed advantages and disadvantages of the extensions available in the literature and proposed a new extension called multidirect. Compared to Lipschitz optimization algorithms, multidirect guarantees global convergence even if Lipschitz constants are not known or do not exist. The introduced notion of potential Pareto optimality reduces the number of hyperintervals selected for subdivision and speeds up the convergence when compared to algorithms using more elementary selection rules. A preliminary numerical investigation confirms the effectiveness of this algorithm. The combination of potential Pareto optimality with different single objective optimization algorithms and further numerical investigations will be subjects of future research. Acknowledgements Thisresearchisrelated tothethematicresearch areaDEMO(Decision Analytics utilizing Causal Models and Multiobjective Optimization, jyu.fi/demo) of the University of Jyväskylä. Funding Open access funding provided by Universitá degli Studi di Padova within the CRUI-CARE Agreement. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, 123