Abraham Gutierrez

Mickaël Buchet
email: buchet@tugraz.at@s.clair

Sylvain Clair
email: sylvain.clair@univ-amu.fr

Persistent homology to quantify the quality of surfacesupported covalent networks

Keywords: On-surface synthesis, Scanning probe microscopy, Persistent homology, 2D networks, Network quality

Covalent networks formed by on-surface synthesis usually suffer from the presence of a large number of defects. We report on a methodology to characterize such two-dimensional networks from their experimental images obtained by scanning probe microscopy. The computation is based on a persistent homology approach and provides a quantitative score indicative of the network homogeneity. We compare our scoring method with results previously obtained using minimal spanning tree analyses and we apply it to some molecular systems appearing in the existing literature.

Introduction

On-surface synthesis is emerging as a very efficient technique to create organic functional surfaces. [1][2][3][4] In a bottom-up approach, well designed molecular precursors are deposited on a surface and activated to produce original chemical compounds or covalently linked networks.

The most successful realization of on-surface synthesis is the creation of graphene nanoribbons (GNR) with perfectly defined atomic structure. [5][6][START_REF] Corso | Bottom-Up Fabrication of Atomically Precise Graphene Nanoribbons[END_REF] While such one-dimensional (1D) approach remains rather easily accessible and produces polymeric wires up to the micrometer range, [START_REF] Zhong | [END_REF][9][10][11][12] the quest for producing two-dimensional (2D) covalent structures, or 2D polymers, [13] is highly challenging. The formation of surface-supported 2D covalent structures can be seen as an extension of graphene growth, whereby the extraordinary properties of this material [14] are improved further by controlling the atomic scale structure and the introduction of heteroatoms. Robust surfaces with exceptional properties are thus expected. [15][16][17] A large variety of surface-supported 2D networks were fabricated in ultrahigh vacuum (UHV) environment [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] or at the liquid-solid interface. [37][38][39][40][41][42][43][44] In the latter case thermodynamic equilibrium conditions can be reached and extended well-ordered networks can be formed. In contrary, in UHV the growth process is irreversible, which prohibits self-healing mechanisms and generates a large number of defects. Indeed, most 2D networks to date are limited to a few tens of nm in size and are poorly ordered. The formation of defects in the networks can have several origins: 1) incomplete reactions, leading to the formation of oligomeric portions with free ends; 2) flexibility of the molecular backbone. Some distortions in the network can arise from slight variations in the bond geometries as compared to the ideal cases. For example in honeycomb networks pentagonal and heptagonal pores, or even tetragonal and octagonal pores are sometimes observed alongside the ideal hexagonal pores [20,21,[45][46][47] ; 3) poor selectivity in the reaction mechanism, i.e. the introduction of different products with poorly adapted configuration. A wide set of tunable parameters is usually available to direct the network growth, such as the reaction temperature, the substrate nature and its crystallographic orientation, or the addition of metal adatoms. [4] However, the complexity level in on-surface synthesis is high and it is still hardly possible to extract any general tendency for the reaction mechanisms. [4,[START_REF] Björk | Kinetic and Thermodynamic Considerations in On-Surface Synthesis[END_REF][START_REF] Björk | [END_REF] Nevertheless, the optimization of the growth process and the quest for networks of ideal or well-controlled quality is of prime importance.

Surface-supported networks are usually characterized by scanning probe microscopy (scanning tunneling microscopy, STM, or atomic force microscopy, AFM) The development of an easy and straightforward quantification tool of these real space images is required.

In fact the assessment of the quality of the networks has been barely addressed and the successful network formation is usually qualified subjectively as of "good" or "poor" quality.

The different defects observed can be assigned from the STM images and counted individually to get an estimation of the overall quality. Only the extension of the networks (domain size) can be properly quantified. The global morphology of the networks can be assessed qualitatively by using Monte Carlo simulations and their comparison with the STM images. [50] A statistical analysis approach based on the minimal spanning tree (MST) method [51] to provide a quantitative and comparative estimation of the quality of boronic acid based covalent networks has been proposed. [45] This approach represents, to our knowledge, the one and only attempt to treat the problem of the quantification of the network quality.

Persistent homology has been recently proposed as an efficient topological method to analyze the structure of materials. For example, it was used for the analysis of pore configuration of granular materials, [52] for qualitative discrimination of amorphous metals [53] or for an analysis of similarity between nanoporous materials. [54] In this work we used an approach based on persistent homology to provide a simple and straightforward tool assessing the network quality in a quantitative way. Persistent diagrams were produced from STM topographic images and analyzed to provide a numerical score. We show that our method is consistent with the previously published MST approach and is robust against the image quality. The method was further used with a few additional results from the literature to highlight its general usability.

Results and discussion

Persistent homology has been used in previous applications in material science and interested readers can refer to Ref. [START_REF] Buchet | Persistent Homology and Materials Informatics[END_REF] as an introduction to the methodology. A persistent diagram was computed from the network represented in the experimental image. Persistent diagrams represent the evolutions of holes (the pores) while the line width inside the network is artificially varied. Each hole in the structure is given a birth time when the enclosing circle first appears and a death time when the hole is completely filled. Each point in the persistent diagram corresponds thus to the pair (birth, death) for a specific hole (see Fig. 1). Details of the method are provided in supplementary information. The code used is available for download. [56] In the persistence diagram the columns at birth time ≤ 0 correspond to the pore size distribution of the image.

Contrary to most applications of persistent homology, we are not concentrating on the holes with the largest lifespan. Instead we want to evaluate how regular the network is. Therefore Intuitively, the PH score gives the proportion of holes similar to the most representative hole. Therefore the higher is the score, the more regular the network is. A score of 1 would be achieved for a perfectly regular network, independently of the form of the network under the conditions that only one type of cell exists and that these cells are convex. The method is robust with respect to boundary effects. Holes formed along the boundary do occur in certain cases but their number is of a different order of magnitude than the total number of holes.

Therefore they only modify marginally the score. For this reason, PH scoring is poorly affected by the domain sizes, the surface coverage, or the molecular concentration. The method is also robust with respect to the image size (see Supplementary Information Fig. S1), and images taken at different scales can be comparatively analyzed. First we focused on the work by Ourdjini et al. [45] that provided a systematic analysis of the quality of the networks obtained from the self-condensation of 1,4-benzenediboronic acid (BDBA) on coinage metal surfaces. The different growth conditions were classified according to the outcome of the MST analysis. In brief, the centers of all pores were connected by straight segments, and the standard deviation of the segment lengths (MST ) is plotted versus the mean segment length (MST m). The results from the two methods (MST and PH) could be directly compared and are reproduced in Fig. 2. We found thus that the PH score, that considers the inner pore sizes, is consistent with the MST result that considers the interpore distances. The highest score obtained with BDBA, on Ag(111) surface, was given a PH score of 0.34.

An improvement of the network quality was proposed in Ref. [46] by using a sequential growth strategy [4] with the modified precursor p-bromobenzene boronic acid (BBBA) molecule on Au(111). The MST analysis is reported in Fig. 2. The higher quality gained with this precursor was confirmed by the PH score of 0.53.

All the above mentioned networks were grown in UHV conditions, in which the selfcondensation reaction is a strongly irreversible process, thus leading to overall rather low quality of the networks. Equilibrium growth conditions were reached in atmospheric conditions on a HOPG surface, while the network was formed by annealing the surface in the presence of water pressure confined in a closed reactor. [38,41] In this case (see Fig. 2), a nearly perfect honeycomb network was formed with the larger 4,4'-biphenyldiboronic acid (BPDA) precursor, [38] and for which we calculated a PH score of 0.88. Comparatively, in similar conditions but without the presence of H2O in the reactor, [38] a lower score of 0.20 is measured (Fig. 2).

Recently the smallest network based on boronic acid chemistry was achieved using the most basic diboronic acid precursor, namely tetrahydroxydiboron (THDB, see Fig. 2). [24] However, several defects are formed in this network, and our PH analysis gave a rather low quality score of 0.18. All the results obtained with boronic acids are summarized in Table 1. In fact a large pore wall length (i.e. a longer diboronic acid precursor) can provide some geometrical flexibility inducing the formation of defective pores with non-ideal shapes (deviating from the hexagonal symmetry), even if the reaction is complete and all the boronic acid groups have reacted. It was thus shown in Ref. [46] that, for networks made from BBBA precursors and comprising biphenyl groups as pore walls, polygonal pore shapes ranging from tetragons to octagons were routinely formed, whereas only pentagons to heptagons were observed using BDBA precursors, for which the pore wall length is limited to one phenyl unit. [35] This observation is however not a general rule as also tetragonal pores were observed with THDB. [24] We can conclude from the comparative analysis of the PH scores for these different networks in UHV that the network quality is not obviously affected by the geometrical flexibility of the covalent links. Networks #1 to #6 are obtained with BDBA precursor, see Table 1 and Ref. [45] for experimental details. For networks obtained with BBBA, BPDA and THDB, see Refs. [46] , [24] , [38] , respectively, for experimental details.

Adapted with permission from [45] Copyright 2011 by the American Physical Society.

Adapted with permission from [46] . Copyright 2012 American Chemical Society. Republished with permission of Royal Society of Chemistry, from [24] and [38] ; permission conveyed through Copyright Clearance Center, Inc.

Precursor / surface Annealing temperature PH score Reference BDBA / Au(111) #6 250°C 0.09 [45] BDBA / Ag(100) #5 25°C 0.13 [45] BDBA / Cu(111) #4 150°C 0.16 [45] BDBA / Au(111) #3 25°C 0.20 [45] BDBA / Ag(100) #2 25°C 0.27 [45] BDBA / Ag(111) #1 25°C 0.34 [45] BBBA / Au(111) 250°C 0.53 [46] THDB / Au(111) 25°C -175°C 0.18 [24] BPDA / HOPG 0.88 [38] equilibrium growth BPDA / HOPG ambient growth without H2O 0.20 [38] Table 1. Summary of PH scores obtained from boronic acid based networks. Adapted with permission from [50] . Copyright 2010 American Chemical Society.

Now we consider the formation of so-called porous graphene networks obtained from the hexaiodo-substituted macrocyclic cyclohexa-m-phenylene (I6-CHP) precursor. [34,50] A comparative study of the network growth on the Cu(111), Ag(111) and Au(111) surfaces was proposed in Ref. [50] . The STM images were compared with Monte Carlo simulations whereby the relevant parameter P was shown to be the ratio between the reactivity (coupling probability) and the diffusivity. A high P value was proposed to be valid for the growth on Cu(111), an intermediate P for Au(111) and a low P for Ag(111) (see Fig. 3). For these three surfaces we found PH scores of 0.22, 0.35 and 0.48, respectively, thus confirming Ag(111) as the most adapted surface for an ideal network growth. Note that for this specific system the molecular precursor is imaged with a donut shape exhibiting an internal pore. The latter is contributing to the final score and in particular to the non-negligible value of 0.22 for the Cu(111) case. Nevertheless, the increase in quality and the extension of the porous graphene network are well reproduced in the PH score evolution. See text for details.

(b) Republished with permission of Royal Society of Chemistry, from [21] ; permission conveyed through Copyright Clearance Center, Inc. (c) Reprinted from [19] , CC BY 4.0. (e,f)

Republished with permission of Royal Society of Chemistry, from [26] ; permission conveyed through Copyright Clearance Center, Inc.

The honeycomb network formed by the precursor 1,3,5-tris(4-bromophenyl)benzene (TBPB, Fig. 4a) has been studied by different groups. On Au(111), a nicely ordered extended network was obtained, [21] for which we found a PH score of 0.53 (see Fig. 4b). Note that here the score is lowered due to the presence of filled pores that are thus wrongly assigned as defects by the calculation. On a Au-enriched Pd(111) surface, after annealing to 400°C a disordered network was obtained, for which we calculated a PH score of 0.19. [19] Further annealing up to 510°C allowed to improve the network quality and the PH score raised to 0.30 (Fig. 4c). The same honeycomb network could be obtained using a sequential coupling strategy [4] with a precursor bearing two different halogen types (4,4''-dichloro-1,1':4',1''-terphenyl, DCTP, Fig. 4d). [26] The low quality network obtained on Au(111) with the help of codeposited Cu adatom catalysts delivered a PH score of 0.16 (Fig. 4e). Similarly, on Cu(111) a score of 0.14 was measured (see Fig. 4f). [22] Copyright 2010. (b) Adapted with permission from [9] . Copyright 2014 American Chemical Society.

Finally, we present two additional systems to show the generality of our method. Most of the covalent networks reported in the literature have a honeycomb-like structure, but an extended network with an ideally square symmetry was reported with the molecule 5,10,15,20tetrakis(4-bromophenyl)porphyrinato}nickel(II) (NiTBrPP) on the Au(111) (see Fig. 5a). [22] Our PH scoring method could be applied in the same way on this network and delivered a score of 0.24. To end with, we show the example of a very low quality network obtained with the molecule 4,4'-diethynyl-1,1':4',1''terphenyl on Ag(111). [9] In this case, in addition to threefold linkages, twofold linkages between the linear precursors are possible with similar occurrence probability, which obviously prevents the formation of an ordered network.

Nevertheless, here the 2D character of the network is preserved (Fig. 5b). Due to the high number of defects, such network can be considered as a network with very low ordering.

Indeed, a PH score of 0.07 was calculated for it.

Conclusion

We propose a new method based on persistent homology to assess the quality of surfacesupported covalent networks from STM images. We show that our method is consistent with the previously published tool based on the minimal spanning tree (MST) approach. [45,46] The PH score provides quantitative insights on various types of networks, an issue that remains of prime importance for developing a rational optimization of the growth conditions in onsurface synthesis.

The scoring for the regularity of the networks relies on using a new approach to persistent homology, potentially discarding high lifespan topological features to concentrate on the density of the diagram. We believe that this approach could be also useful for the quantification of a wide range of other systems organized in networks and characterized by 2D imaging, such as block copolymers, [START_REF] Segalman | [END_REF]58] but also 2D soap froth, [59] biological cells, [60] Bénard-Marangoni convection cells, [61] or even geological structures. [62] Also, there is no theoretical objection to use it for higher dimensional structures. For example, in a given a set of slices representing a 3D system it could be possible to analyze the regularity of the structure using the cavities instead of the holes.

  we are interested in how concentrated the persistence diagram is. A score is computed for each point by considering the number of similar points within a certain window, divided by the total number of points in the diagram. The score of the whole diagram (called in the following the PH score) is the maximal such score that can be obtained across the points in the diagram.

Fig. 1 .

 1 Fig. 1. (a) Representation of a simulated network with honeycomb structure for which defects are randomly introduced. (b) Corresponding persistent diagram. The multiplicity corresponds to the number of occurrences of each point. A score of 0.75 was calculated for the network.

Fig. 2 .

 2 Fig. 2. Surface-supported covalent network formed by self-condensation of boronic acids. The computed PH scores are compared with results obtained by minimal spanning tree analyses.

Fig. 3 .

 3 Fig. 3. Evolution of the PH score with different substrates for the growth of porous graphene.

Fig. 4 .

 4 Fig. 4. PH scores obtained for various honeycomb networks formed by Ullmann coupling.

Fig. 5 .

 5 Fig. 5. (a) PH score of a network with intrinsic square symmetry obtained by Ullmann coupling. (b) PH score of a highly disordered network obtained by homocoupling of terminal alkynes. (a) Reprinted by permission from Springer Nature: Nano Research, Ref.[22]

Acknowledgements

Michael Kerber and Yasuaki Hiroaka are acknowledged for motivating discussions. A. G. acknowledges financial support from the Austrian Science Fund project FWF P29355-N35.

Financial support from Centre National de la Recherche Scientifique (CNRS), Aix-Marseille University and Agence Nationale de la Recherche (ANR grant N° ANR-17-CE08-0010 "DUALITY") is acknowledged.


Nicolas Berveglieri
email: nicolas.berveglieri@univ-lille.fr

Bilel Derbel
email: bilel.derbel@univ-lille.fr

Arnaud Liefooghe
email: arnaud.liefooghe@univ-lille.fr

Hernan Aguirre

Kiyoshi Tanaka
email: ktanaka@shinshu-u.ac.jp

Hernán Aguirre

Surrogate-assisted multiobjective optimization based on decomposition

Keywords: Theory of computation → Gaussian processes, Stochastic control and optimization, Algorithm Multiobjective optimization, surrogates, benchmarking

come

INTRODUCTION AND BACKGROUND

General Context. Multiobjective optimization problems (MOPs) [START_REF] Tamaki | Multi-objective optimization by genetic algorithms: a review[END_REF][START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative case study and the strength pareto approach[END_REF] aims at simultaneously optimizing two or more objectives. Given that these objectives are typically in conlict, there exist a whole set of optimal solutions that provide diferent quality tradeofs. The challenging task of computing a high-quality approximation set can be tackled using multiobjective optimization evolutionary algorithms (MOEAs). MOEAs were proved to be well-suited for blackbox problems, for which a mathematical formulation of the objective functions might not be available. In this paper, we are motivated by the challenges underlying the high cost of computing the objective values, e.g., when running a heavy simulation process [START_REF] Branke | Multiobjective optimization, interactive and evolutionary approaches[END_REF]. More speciically, we are interested in an expensive optimization scenario [START_REF] Tenne | Computational Intelligence in Expensive Optimization Problems[END_REF], where the cost of computing the objective value(s) impose drastic restrictions on the number function evaluations. Depending on the CPU time that can be aforded, an algorithm designer must accommodate to a relatively short budget, typically few hundreds to few thousands evaluations. As a consequence, using standard evolutionary operators, which are supposed to serve the search in a non-expensive setting, become problematic because of the diiculty of inely controlling their stochastic behavior. In this context, adapting conventional MOEAs so that they can operate in a combined manner with well-established machine-learning techniques have attracted much attention in recent years [START_REF] Deb | A taxonomy for metamodeling frameworks for evolutionary multiobjective optimization[END_REF][START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF].

Related work overview. On the one hand, a general idea that was extensively explored so far is to build one (or multiple) surrogate model(s) from which some information about the (unknown blackbox) objective functions can be extracted and injected to the search process in order to help reaching the most interesting regions more eiciently [START_REF] Loshchilov | Surrogate-Assisted Evolutionary Algorithms[END_REF]. A broad range of investigations can be found on using surrogates to assist MOEAs, such as those based on support vector regression [START_REF] Pavelski | Elmoea/d-de: Extreme learning surrogate models in multi-objective optimization based on decomposition and diferential evolution[END_REF], regression trees [START_REF] Sagawa | Learning variable importance to guide recombination[END_REF], radial basis functions [START_REF] Saúl | MOEA/D assisted by rbf networks for expensive multi-objective optimization problems[END_REF], or Gaussian processes [START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with gaussian process model[END_REF]. Depending on the nature of the considered surrogates, several design options might be adopted. This ranges from standard approaches where cheap surrogates are used in lieu of the expensive objective functions, to those using the surrogates for pre-screening solutions generated by the evolution engine. Other options from single-objective Bayesian global optimization and reinforcement learning consist in a criterion to inely tune the exploitation-exploration trade-of when sampling new points. A number of merit functions and inill criteria were already investigated and combined with MOEAs for selecting the most beneicial solutions to be evaluated [START_REF] Jones | Eicient global optimization of expensive black-box functions[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF].

On the other hand, despite the substantial progress achieved in the last decade, there are relatively few systematic investigations aiming at providing a better understanding of what makes a particular approach distinguishable from others, and under which optimization scenario, in terms of budget and problem diiculty, it is recommended. We can cite [START_REF] Deb | A taxonomy for metamodeling frameworks for evolutionary multiobjective optimization[END_REF] which provides a coarse-grained taxonomy of surrogate-assisted MOEAs with a particular focus on constrained MOPs. In [START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF], one can also ind a complementary taxonomy, and a discussion on extending existing approaches to support the production of a batch of multiple points that can be evaluated in parallel. Besides considering a simple parallel extension of ParEGO [START_REF] Knowles | Parego: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF], using a set of weight vectors, the most advanced MOEAs studied therein fall into the class of dominance-based and indicator-based algorithms such as SMS-EGO [START_REF] Ponweiser | Multiobjective optimization on a limited budget of evaluations using model-assisted s-metric selection[END_REF]. This taxonomy came with an interesting empirical analysis of the diferent considered extensions with a focus on the potential gain in terms of parallelism. Nonetheless, a class of state-of-the-art EMO approaches was only considered at a small extent, namely, surrogate-assisted variants of the so-called MOEA based on Decomposition (MOEAD) [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF].

The MOEAD (multiobjective evolutionary algorithm based on decomposition) framework belongs to the class of aggregationbased approaches. The original MOP is transformed into a number of (single-objective) sub-problems, deined using a scalarizing function, that are solved cooperatively. Due to its eiciency and lexibility in integrating existing single-objective optimizers into the multiobjective setting, MOEAD has become one of the most studied approaches. Nonetheless, when dealing with an expensive setting, we are not aware of any systematic analysis of what makes decomposition beneicial and how to efectively integrate surrogates in the original cooperative solving process of MOEAD. Interestingly, the most recent surrogate-assisted approaches are based on leveraging concepts from multiobjective decomposition [START_REF] Chugh | A surrogate-assisted reference vector guided evolutionary algorithm for computationally expensive many-objective optimization[END_REF][START_REF] Hussein | A generative kriging surrogate model for constrained and unconstrained multi-objective optimization[END_REF][START_REF] Pavelski | Elmoea/d-de: Extreme learning surrogate models in multi-objective optimization based on decomposition and diferential evolution[END_REF][START_REF] Saúl | MOEA/D assisted by rbf networks for expensive multi-objective optimization problems[END_REF][START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with gaussian process model[END_REF]. This is precisely where the contribution of this paper lies.

Methodology and contributions overview. This paper can be viewed as the continuation of previous works on understanding and delimiting the core concepts of surrogate-assisted MOEAs [START_REF] Deb | A taxonomy for metamodeling frameworks for evolutionary multiobjective optimization[END_REF][START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF], with a particular focus on decomposition techniques. We thereby propose to discuss a ine-grained classiication of existing approaches in light of the taxonomy provided in [START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF]. This is in an attempt to grasp the speciicity of existing decomposition-based MOEAs, and to highlight the common algorithmic components that can be integrated interchangeably in diferent approaches. We also include simple archetype algorithms serving as baseline and which are barely considered in previous works. In fact, besides reviewing a number of distinguishable surrogate-assisted MOEAs within the framework of existing taxonomies, our goal is to shed the light on the strengths and weaknesses of existing approaches through an extensive benchmarking efort. Indeed, unfortunately, despite a skillful design, existing algorithms are often validated on diferent benchmark functions and under diferent, sometimes incomparable, settings. In our work, we aim at providing a comprehensive feedback on the relative performance of advanced approaches from the literature, especially when compared to a priori simple design choices. We rely on a number of problems taken from the bi-objective black-box optimization benchmarking test suite (bbob-biobj), and providing a representative sample of single-objective function combinations with well-known facets of diiculty [START_REF] Mersmann | Coco: A platform for comparing continuous optimizers in a black-box setting[END_REF]. Our empirical indings are to be considered as a step towards providing a uniied view of surrogate-assisted MOEAs based on decomposition, and their relative performance under diferent scenarios. In particular, we study the impact of clustering when training the surrogate models, both on running time and on quality, as well as the impact of problem characteristics on convergence and running time.

Outline. In Section 2, we provide a classiication of some distinguishable surrogate-assisted MOEAs. In Section 3, we report our experimental indings. In Section 4, we conclude the paper.

SURROGATE-ASSISTED MOEA BASED ON DECOMPOSITION (S-MOEAD)

In this section, we provide a ine-grained classiication of existing surrogate-assisted MOEAs based on decomposition (S-MOEAD) in light of the taxonomy introduced in [START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF]. Our goal is to capture the design choices to be made when combining surrogate models with decomposition, hence eventually leveraging existing frameworks. We shall then provide a more focused description of some representative and state-of-art S-MOEAD algorithms from the literature, illustrating how the ine-grained classiication captures them. The discussed algorithms will also serve to conduct our benchmarking analysis, which is the core contribution of the paper. Before going in a more detailed discussion, let us irst deine the multiobjective decomposition framework in the context of our work.

Decomposition-based MOEA in a Nutshell

Solving a MOP using decomposition implies to deine a number of a priori smaller single-(or multiobjective) sub-problems. In our work, we assume a sub-problem is simply deined by aggregating the original objectives using a scalarizing function. We consider an objective function vector F : I R d → I R m to be minimized, such that I R d is the (continuous) variable space and I R m is the objective space. For two solutions x, x ′ ∈ I R d , x is dominated by x ′ if, for all i ∈ {1, . . . , m}, f i (x) ⩽ f i (x ′ ), and there is a j ∈ {1, . . . , m} such that

f j (x) < f j (x ′ ). A solution x ⋆ ∈ I R d is Pareto optimal if it there does not exist x ∈ I R d s.t. x ⋆ is dominated by x.
The set of all Pareto optimal solutions is the Pareto set; its mapping in the objective space is the Pareto front. Let w ∈ I R m be a weight vector, and д(•|w) be a scalarizing function that, given a solution, transforms its objective vector F (x) into a scalar real value. For instance, one might want to use the Chebyshev function, deined by

д(x |w) = max i {w i •| f i (x)-z * i |}
, where z * is a reference point. Other scalarizing functions can be found in the MOEA literature [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF].

To identify a good Pareto set approximation, a standard decomposition approach considers to solve a number of sub-problems, deined using diferent weight vectors for coniguring the scalarizing function. The output is hence made up with diferent solutions, corresponding to the deined sub-problems. Sub-problems may be solved independently [START_REF] Hughes | Multiple single objective Pareto sampling[END_REF] or cooperatively [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF]. Besides, diferent evolutionary mechanisms can be adopted to evolve the solution set [START_REF] Trivedi | A survey of multiobjective evolutionary algorithms based on decomposition[END_REF]. In an expensive setting, and looking at how surrogates can be combined with decomposition, most approaches share the use of a scalarizing function. However, diferent components and strategies are adopted to ind good solutions w.r.t. the scalarized sub-problems in a restricted amount of evaluations. This makes it relatively diicult to diferentiate between existing approaches and to identify their key components at irst sight. In the next section, we describe a simple ine-grained classiication of S-MOEAD in an attempt to capture their algorithmic gestalt.

A Fine-grained Classiication of S-MOEAD

2.2.1 General outlook. The considered classiication is summarized in Fig. 1. It is closely related to the taxonomy provided in [START_REF] Deb | A taxonomy for metamodeling frameworks for evolutionary multiobjective optimization[END_REF], and is to be viewed as a ine-grained reinement targeting a better understanding of S-MOEAD. In fact, most surrogate-assisted MOEAs are following seemingly the same general template with basically four a priori interdependent layers, in addition to the initialization step. Initialization is particularly important in an expensive setting, since every single evaluation counts. Like all surrogate-assisted approaches, an S-MOEAD algorithm has to specify the set of initial points. This is not only of crucial importance for the search process, but also for the learning process, as the population serves as a pool from which the training set is constructed. Most algorithms use a Latin Hypercube Sampling (LHS), bare of any problem-speciic heuristic. The general idea captured by the other layers is then as follows. Given a set of points from which the actual objective values are known, a pre-processing mechanism is used to prepare the data for training. Having built the training set, one or several models are itted for a well-deined model response. These models are intended to help the search process by guiding the sampling and evaluation of new solutions in the subsequent layers. It is important to understand that the generation of new candidate solutions and the selection of solutions for evaluation are two diferent aspects that should be distinguished. In fact, a critically-important step is to generate a set of candidate solutions explicitly based on the surrogate models, or on the contrary implicitly by allowing the surrogate models to interfere with the conventional evolutionary search process. This is inally followed with a selection step that decides which subset of solutions is to be evaluated using the actual expensive objectives. We argue that the speciication of these layers and the setting of their components allows one to diferentiate a particular approach form another. It is worth noticing that the initial taxonomy provided in [START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF] does not include the training set construction phase, which will be shown to be a key ingredient in diferent S-MOEAD approaches. In the following, we further discuss the diferent layers in a more ine-grained manner, before turning to a more focused discussion of existing techniques.

Defining the model training set.

This phase can take diverse forms, as almost each existing algorithm adopts a diferent approach. It is crucially important since the input data can highly bias the accuracy of the constructed model and/or inluence the interpretation of any information extracted. A straightforward strategy, not necessarily speciic to S-MOEAD, is simply to include all solutions evaluated so far as input for the model training/itting and for the subsequent model-speciic treatments. However, this might lead to the construction of a global surrogate at a high computational cost, eventually dominate the evaluation cost itself. A common alternative is to use only recently-evaluated points, that can be selected, e.g., over a predeined time window, or also based on the quality of the observed objective or scalar values. Another alternative is to construct the training set using a clustering procedure. This is a natural outcome in S-MOEAD approaches, since dealing with diferent sub-problems suggests to focus on speciic target search regions directed by the value of the weight vector, and hence to clusterize the data accordingly. Interestingly, this could imply the use of a structured ensemble of surrogates used in the other layers to assist the optimization process. As such, we can ind a number of data clustering techniques, ranging from the ones operating in the variable space, to those centered on the objective values of the points sampled so far, or even on the aggregated values deined with respect to the scalarizing function and/or its weight vector. From a pure machine learning perspective, choosing/clustering the evaluated solutions is both related to the degree of locality expressed by the trained model and to the possibility of obtaining diferent trade-ofs between the model quality and the computational cost of the itting phase [START_REF] Snelson | Local and global sparse gaussian process approximations[END_REF].

Model fiting and response.

At this stage, one is concerned with the output prediction provided by the surrogate model, i.e. the response variable. A number of variants can be considered, depending on the type of training data and on the surrogate model by itself. The most straightforward approach is to it the objective values, meaning that one model is itted w.r.t. each objective. The scalarizing function can then be used to aggregate the predicted values. The other common approach is to it the aggregated values, and an ensemble of models is hence constructed w.r.t. the weight vectors of the scalarizing functions. This might then be tightly coupled with a speciic choice of the training set, eventually leading to a complex design. Besides, since diferent regression models and diferent hyper-parameters can be considered, diferent models can actually be trained concurrently and their responses exploited afterwards.

2.2.4 Candidate solution(s) generation. This phase consists in computing a set of solutions that are expected to be the most beneicial at the current state of the search process. This is critically important, since it speciically informs about the most promising regions that should be efectively searched. We identify three main aspects to this phase. First, one can target candidate solutions optimizing explicitly the model response, i.e. the objective values or their aggregation as predicted by the surrogate model(s). In this case, a standard evolutionary algorithm (EA) can for instance be carried out, while using the model response as a substitute for the real objective or scalarizing functions.Moreover, we diferentiate between approaches running the EA for multiple iterations, e.g., run a whole MOEAD up to some predeined number of generations, or those running a single EA iteration before resuming, e.g., generate multiple ofspring using standard variation operators and then choose a subset as candidate solutions. The former approach is usually termed as substitute approach [START_REF] Pavelski | Elmoea/d-de: Extreme learning surrogate models in multi-objective optimization based on decomposition and diferential evolution[END_REF], whether the latter is usually termed as iltering or pre-screening approach [START_REF] Loshchilov | Surrogate-Assisted Evolutionary Algorithms[END_REF]. Second, one can target candidate solutions optimizing other criteria that are intrinsically linked with the surrogate model. This is typically the case of Gaussian processes, where a number of inill criteria are well established, e.g., Expected Improvement (EI), Probability of Improvement (PI), or Lower Conidence Bound (LCB) [START_REF] Auer | Using conidence bounds for exploitation-exploration trade-ofs[END_REF][START_REF] Hawe | Probability of improvement methods for constrained multi-objective optimization[END_REF][START_REF] Kleijnen | Expected improvement in eicient global optimization through bootstrapped kriging[END_REF]. However, since we are dealing with a MOP, one could for instance search for a single solution optimizing the inill criteria with respect to one or multiple weight vectors, or run a standard MOEA where the objectives are deined as the inill criteria corresponding to each original objective. Finally, these diferent design choices can be guided by whether a single candidate solution is to be evaluated, or if a batch of solutions can be evaluated in parallel. Further considerations on this aspect are discussed in [START_REF] Deb | A taxonomy for metamodeling frameworks for evolutionary multiobjective optimization[END_REF], not necessarily within the context of decomposition.

Selection for evaluation.

When multiple candidate solutions are generated, a (multiobjective) subset selection mechanism is to be designed in order to prepare a batch of a desired size on which the true objective values is to be supplied. In the context of S-MOEAD, this can be done on the basis of the scalar values retrieved or computed from the model response, or based on the inill criteria, in addition to other options available from the literature. On the contrary, when a single candidate solution is generated while a batch of solutions is expected, one has to iterate over the previous phase, eventually injecting some knowledge about the solutions that are already contained in the batch. Once the batch is ready, the selected solutions are evaluated, and the whole process is repeated until some termination condition is satisied. A maximum number of evaluations is typically adopted as a stopping condition.

Considered S-MOEAD Algorithms

Based on the previous discussion, we are able to classify a number of existing S-MOEAD techniques from the literature, and to diferentiate them in a more systematic manner, as depicted in Table 1. Apart from the four irst approaches, which follow a self-designed simple iltering archetype framework, the remaining one are based on more sophisticated components. Notice that all these approaches can be correctly classiied according to the diferent layers provided in Fig. 1. As one can see, the considered approaches provide a representative combinations of the possible design choices. In the following, we shall describe them in more details and clarify any important hidden technical detail, since we also aim at conducting a comparative analysis of their behaviors in the next section.

S-MOEAD based on filtering.

The approaches based on iltering are conceptually among the most basic S-MOEADs; see, e.g., [START_REF] Siang | Machine learning enhanced multi-objective evolutionary algorithm based on decomposition[END_REF]. Instead of generating one single solution by means of variation operators, a pool of λ > 1 solutions are generated. This pool constitutes the set of candidate solutions from which a batch has to be selected for evaluation. Once the evaluation is carried out, the training set is updated, the model is re-trained, and the next iteration of the original MOEAD can be executed. It should be clear that, even for this very simple archetype, diferent variants could be designed on the basis of the previously-discussed taxonomy. Interestingly, this archetype allows us to freely investigate the impact of setting the components from other layers, such as the impact of clustering or the response variable of the surrogate. We thereby consider four simple ilter-based variants, based on the conventional MOEAD, as depicted in Table 1. The two irst variants employ one surrogate model per objective, whereas the third and fourth ones train one surrogate model for each sub-problem deined by the weight vectors. The irst and fourth variants use the whole set of solutions evaluated so far for training, whereas the second and third ones use two diferent clustering techniques for constructing the training set. The second variant performs a fuzzy clustering w.r.t. the position of solutions in the variable space, and train one model per cluster and per objective. In MOEAD-ilter2, the estimation is based on the model linked to the cluster with the closest center to the considered solution. On the contrary, MOEAD-ilter3 performs a clustering w.r.t. the aggregated values, each cluster being used for a given weight vector. Taken separately, these two clustering approaches are not our own proposal, since they are inspired by key components used in more advanced (state-of-the-art) non-iltering based S-MOEAD [START_REF] Hussein | A generative kriging surrogate model for constrained and unconstrained multi-objective optimization[END_REF][START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with gaussian process model[END_REF]. The detailed description of these clustering techniques will be presented in more details below.

S-MOEAD based on substitution.

A second class of S-MOEAD archetypes is based on using the surrogate as a substitute, by temporarily relying on its response for searching improving solutions.

In [START_REF] Pavelski | Elmoea/d-de: Extreme learning surrogate models in multi-objective optimization based on decomposition and diferential evolution[END_REF], a conventional MOEAD optimizing the model response is executed repeatedly for a number of generations before a batch of evolved solutions is evaluated using the expensive objective functions. However, at each MOEAD iteration, i.e., when looping over weight vectors, a iltering approach selects the solution considered for the replacement phase among a number of ofspring generated by variation. Another conceptually more sophisticated variant based on substitution is described in [START_REF] Saúl | MOEA/D assisted by rbf networks for expensive multi-objective optimization problems[END_REF]. The so-called MOEAD-RBF constructs the training set using the best subset of evaluated solutions w.r.t. the weight vectors. Next, not one, but three radial basis function (RBF) surrogate models per objective are trained using diferent kernels. After training, the objective values from solutions evaluated so far are compared against the corresponding predicted values given by each of the so-constructed RBF model. The observed error is used to rank the accuracy of each model; see [START_REF] Saúl | MOEA/D assisted by rbf networks for expensive multi-objective optimization problems[END_REF] for details. Then, a standard MOEAD is run for a number of generations using the RBF models as a substitute, but where the predicted values are weighted by a linear combination of the previously-computed model weights. The approximation set returned by MOEAD is then used as the set of candidate solutions from which a batch of solutions is to be selected. This is performed by considering a number of weight vectors from which a desired fraction is considered in a round-robin fashion. The candidate solutions having the best model-weighted aggregated values for the selected weight vectors is considered for the expensive evaluation.

S-MOEAD based on

Gaussian processes. Many surrogateassisted MOEAs are based on the extension and adaptation of the popular EGO (Expensive Global Optimization) procedure from single-objective optimization [START_REF] Jones | Eicient global optimization of expensive black-box functions[END_REF]. Existing S-MOEADs stand for no exception. EGO builds a Gaussian Process (GP) model of a function, where typically each point of the search space is viewed as a random variable that is assumed to follow a normal distribution. The generation of promising candidate solutions is fully guided by the knowledge that can be extracted from the model response surface. Diferent merit functions and inill criteria can be designed for scalar (single-objective) functions, such that the widely acknowledged Expected Improvement (EI) [START_REF] Jones | Eicient global optimization of expensive black-box functions[END_REF]. In the following, we review three reference EGO-like S-MOEADs and discuss their speciic design components, in light of our classiication. ParEGO [START_REF] Knowles | Parego: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF] is among the very irst extensions of EGO for MOPs. It is based on the aggregation of the objectives using an augmented Chebyshev scalarizing function. At each iteration, one weigh vector is thrown randomly, and a single-objective EGO is considered for the so-deined scalar sub-problem. A speciic design choice is to train the corresponding model using only the n (user deined) best-ranked solutions w.r.t. their aggregation values. It should be clear that generating weight vectors at random may be problematic, which in our opinion is the major reason motivating several subsequent S-MOEAD algorithms.

A sophisticated extension of EGO in combination with MOEAD can be found in [START_REF] Zhang | Expensive multiobjective optimization by MOEA/D with gaussian process model[END_REF]. The so-called MOEAD-EGO uses a fuzzy clustering procedure to build multiple clusters of solutions evaluated so far. The size of each cluster is set to a user-deined value k. Each cluster is used to build one surrogate model w.r.t. to each objective, i.e., (k • m) models are built at each iteration. A set of candidate solutions are then generated by means of an MOEAD execution to optimize simultaneously the EIs w.r.t. objective. However, it is shown that, when considering a single target weighted scalarizing function, a local surrogate model can be inferred analytically from the initial (non-aggregated) models. The aggregated EIs can then be inferred as well, which makes the optimization process more efective. More speciically, once the aggregated EI values of a candidate solution is required by the MOEAD inner optimizer, the surrogate model corresponding to the cluster with the closest center in the variable space is considered.

More recently, a seemingly diferent and high-performing EGOlike S-MOEAD was proposed [START_REF] Hussein | A generative kriging surrogate model for constrained and unconstrained multi-objective optimization[END_REF]. Contrary to ParEGO, M-EGO considers a pre-deined set of weight vectors, corresponding to diferent sub-problems. It iterates over the weight vectors following a well-deined sequence. At a given iteration, considering a reference weight vector, the training set is chosen as the n closest points to that weight vector, using the projected euclidean distance of all solutions evaluated so far. A local surrogate for the target search direction is then built accordingly, and the optimization of the inill (EI) criteria is performed by means of a scalar EA.

EXPERIMENTAL ANALYSIS 3.1 Experimental Setup and Methodology

Apart from the 8 surrogate-assisted approaches listed in Table 1, we consider 3 surrogate-less algorithms, namely a standard MOEAD, a random search (RS), and a simple Latin Hypercube Sampling (LHS). In all MOEAs, the population size is set to 50, and is initialized by means of LHS. The scalarizing function is Chebyshev, whose reference point is updated with the best objective values seen so far. For iltering algorithms, we use the same parameters as in the MOEAD. At each iltering step, 8 ofsprings are generated the one with the best aggregated predicted value for the current weight is evaluated. We measure the performance of the competing algorithms under diferent scenarios in terms of budget, measured as In terms of problems, we consider a number of numerical functions extracted from the bi-objective black-box optimization benchmarking test suite (bbob-biobj) [START_REF] Mersmann | Coco: A platform for comparing continuous optimizers in a black-box setting[END_REF]. More particularly, we selected 12 functions with diferent properties in terms of separability, conditioning, multimodality, and global structure. They are listed in Table 2. In order to measure the ability of the considered approaches to cope with the problem dimensionality, we experimented functions with d ∈ {2, 3, 5, 10} decision variables. We perform 10 independent runs for each problem and algorithm (5 280 experiment).

The algorithms are compared in terms of the hypervolume relative deviation to the best-found approximation set for the considered function. The hypervolume [START_REF] Auger | Theory of the hypervolume indicator: Optimal-distributions and the choice of the reference point[END_REF] measures the area covered by an approximation set and enclosed by a reference point. Due to the target of our experimental analysis, we consider two diferent settings for the reference point: a "global" reference point that maps to the worst objective value observed on a given instance once the algorithm starts at 50 evaluations, and a "local" reference point that maps to the worst objective value observed on all Pareto set approximations found for a given budget. that for Table 3, the rank obtained are the same for the 2 hypervolume considered.

Overall Approximation Quality

We irst summarize our indings using two scenarios in terms of the global budget: a particularly tight setting of 500 evaluations, and a moderate setting of 2 500 evaluations. Table 3 reports the rank of every algorithm for each pair of functions and problem dimension. Given that 10 algorithms are compared, ranks range from 0 to 9, a lower value being better. The average rank of a given algorithm over all problems and dimensions is reported in the last row. A local hypervolume reference point is used, but we could check that this remains consistent when using a global reference point.

The strategies based on random search (RS) and latin hypercube sampling (LHS) never obtain a better rank than any other approach. Regarding the conventional (surrogate-less) MOEAD, it is almost always outperformed by the other S-MOEAD for a tight budget. However, when the budget is larger, its performance increases, especially when problem dimension is small. It obtains better rank than MOEAD-ilter1, M-EGO and ParEGO when the number of evaluations is 2 500. When comparing the diferent S-MOEAD algorithms, the EGO approaches outperform iltering for a low budget, whereas the opposite holds for a larger budget. Hence, EGO approaches perform better at the early stages of the search process, while iltering approaches seem to slightly more time to converge. More about the convergence proile will be discussed later. Therefore, depending on the budget that can be aforded for the optimization process, one class of algorithms shall be preferred over the other. This highlights the critical importance of the optimization scenario and the careful considerations that must be taken into account while considering the expensive nature of the problem to be solved. It is also interesting to notice that the relative performance of EGO approaches is particularly good for large dimensional problems. For instance, ParEGO is never outperformed by any other approach apart from MOEAD-RBF for d = 10, which is still to be conirmed with even higher budgets. At last, we can highlight the superiority of MOEAD-RBF for both budget scenarios. Indeed, it signiicantly outperforms all other approaches on most instances, with an average rank of 0.4 and 0.2 respectively for 500 and 2 500 evaluations, and is never outperformed by more than two algorithms.

Impact of Clustering

In this section, we investigate the impact of clustering on the performance of S-MOEAD approaches. We recall that clustering enables to select solutions for constructing the training set. An obvious reason for using clustering is to reduce the computational overhead induced by model training, i.e., the larger the training set, the longer the training phase [START_REF] Knowles | Parego: a hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems[END_REF]. On the other side, it is often argued that a larger training set improves the overall model accuracy. However, it remains unclear how clustering actually impacts approximation quality which motivates our further investigations.

Let us focus on iltering approaches, by comparing irst MOEAD-ilter1, where no clustering is performed and where the models learn the objective functions, to MOEAD-ilter2, also learning the objective functions but performing a fuzzy clustering. Similarly, we compare MOEAD-ilter4, where no clustering is performed and where the models are learning scalarized values, to MOEAD-ilter3, where the model also learns from scalarized values, but where a weight-centered clustering is performed. Bear in mind, that beside the clustering step, MEOAD-ilter1 (resp. MOEAD-ilter4) is identical to MOEAD-ilter2 (resp. MOEAD-ilter3). We should also note that the training set for MOEAD-ilter1/4 is made of the entire set of solutions evaluated so far, whereas the cardinality of the training set in MOEAD-ilter2/3 is bounded by 100.

As reported in Table 3, the relative rank of the four algorithms is almost the same for 500 evaluations. We attribute this to the fact that the diference in the training set size is still relatively small at the early stages of the search process. However, for a larger budget of 2 500 evaluations, the diference starts to show, and MOEAD-ilter1 (resp. MOEAD-ilter4) is clearly outperformed by MOEAD-ilter2 (resp. MOEAD-ilter3). In both cases, the algorithms using clustering manage to outperform their respective counterpart. This means that not only clustering allows for a faster computational time, as we will emphasize later, but it also helps to improve solution quality.

Convergence Proile

Convergence is a critical issue, especially in regards of the expensive setting we are interested in. Fig. 2 reports the evolution of the hypervolume relative deviation as the number of evaluations grows (for dimension d = 5). A close look at the relative convergence proile allows us to emphasize how much diferent are iltering approaches from EGO approaches. While the former seem to convergence almost log-log linearly in many cases, the latter reach a particularly good approximation quality just after about 250-500 evaluations, before getting stuck to a convergence point from which only minor improvements, if any, can be observed. Interestingly, the convergence proile of MOEAD-ilter3 and MOEAD-ilter4 seem to follow the same trend than the surrogate-less MOEAD, but with better hypervolume values all the way through. Ultimately, iltering approaches typically end up outperforming EGO approaches, after about 2 000 evaluations. Regarding MOEAD-RBF, it seems to conciliate the best of both worlds. In a similar way than EGO approaches, it converges particularly quickly at the early stages. However, it does not get stuck, and is able to continue improving at a slower convergence rate. In fact, its convergence proile dominates all others for all problems whose irst objective is a Sphere function (f 1 -⋆), and for the most diicult considered problem (f 20 -f 20 ).

Although we have no certainty as to why EGO approaches tend to get stuck at an early stage of the search process, a possible explanation might be the use of the expected improvement (EI) as inill criterion. Indeed, EI is used to discriminate solutions, with the aim of guiding the search process to the most interesting areas while enabling a good balance between solution quality and model reinement. Fig. 3 shows the EI value of the candidate solution f8-f14 f8-f15 f8-f20 f14-f14 f15-f15 f20-f20 f1-f1 f1-f8 f1-f14 f1-f15 f1-f20 f8-f8 #eval hypervolume relative deviation selected for evaluation at each step of a typical ParEGO run on the Sphere-Sphere problem for d = 2. Note that a similar trend is observed for other EGO approaches. We remind that EI values are positive, and the solution with the largest EI value is to be selected.

• • • • • • • MOEAD MOEAD-filter3 MOEAD-filter4 MOEAD-RBF MOEAD-EGO M-EGO ParEGO
Besides, once the model indicates that no improvement is expected, EI makes no diference between a descent and a poor solution, they will both have an EI value of 0. At the initial steps, the maximum EI value increases to a value close to 800 during the irst 350 iterations, before abruptly dropping to values close to 0, apart from small peaks observed at isolated iterations. This means that, after some time, a lot of solutions actually have an EI of 0. By comparing this trend to the hypervolume relative deviation obtained over time, also reported in the igure, we clearly see the impact of EI values on the search progress. This seems to be a law of the EI inill criterion. As noted in [START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF], the performance of ParEGO could be improved using other inill criteria. This suggests that the EI measure can indeed be the source of the observed issue for EGO approaches. Further research on alternative inill criterion or some work-around for the observed issue, might lead to an improvement of EGO approaches.

OUTLOOK AND OPEN ISSUES

In this paper, we investigated surrogate-assisted MOEAs based on decomposition for expensive optimization. Based on a reined classiication, we were able to instantiate existing approaches as well as new alternatives under a common framework, then highlighting their key design choices and components in a more systematic manner. A comparative analysis was conducted on a comprehensive benchmark of bi-objective black-box numerical optimization problems from the bbob-biobj test suite. Our main indings are as follows.

Firstly, it appears clearly that carefully selecting solutions for training the model is of high importance, not only to reduce the amount of data for accelerating the training phase, but also to improve solution quality. Secondly, training multiple models implies to have an ensemble of surrogates. Combined with the design where the training set is clustered around weight vectors, surrogates are then likely to specialize to speciic regions of the objective space, then matching the rationale of decomposition. As such, we argue that a more systematic investigation of this design choice would allow for explicit ways of designing and training specialized local surrogates. This is perfectly in line with the overall good performance of MOEAD-RBF, which also has the particularity of producing a whole batch of solutions to be evaluated. This is certainly of high importance in terms of computational complexity: not only this reduces the number of computationally-demanding model training tasks, this also enables the evaluation of multiple solutions at once in a parallel computing environment. The selection of the batch of solutions to be evaluated by the expensive objective functions shall then carefully cope with the multiobjective nature of the problem at hand by targeting diverse regions of the Pareto front. At last, the behavior of approaches based on Gaussian processes makes them very attractive for an extremely expensive optimization scenario, where the number of evaluations is reduced to the minimum. Notice that basic iltering approaches can be preferred when the application context can accommodate with a restricted but relatively moderate budget. Besides, the use of EI to select solutions seems to quickly inhibit the search progress, so that alternative inill criteria combined with adaptive and hybrid design choices of candidate solutions generation seems to be a promising future research line.

Figure 1 :

 1 Figure1: A reined taxonomy inspired from[START_REF] Horn | Model-based multi-objective optimization: Taxonomy, multi-point proposal, toolbox and benchmark[END_REF] and targeting a ine-grained classiication of S-MOEAD. The text in colour correspond to algorithm components tailored to decomposition-based MOEAs.

Figure 2 :

 2 Figure 2: Convergence proile of the hypervolume relative deviation under a global reference point with respect to the number of expensive evaluations performed, for all problems of dimension d = 5. Notice the log scale on both axes.

Figure 3 :

 3 Figure 3: EI values of evaluated solutions and convergence plot for a typical run of ParEGO on f 1 -f 1 with d = 2.

Table 1 :

 1 S-MOEAD techniques as instances of the general framework.

	algorithm	II ś model training set	III ś model itting and response	IV ś candidate solution(s)	V ś selec. for evaluation
	MOEAD-ilter1	full dataset (1)	1 model per objective (1.1)	1 MOEAD iteration (2.1.1)	single, w.r.t. cur weight (1.2.1)
	MOEAD-ilter2	fuzzy clustering [27] (3.1)	1 model per cluster and objective (1.1) 1 MOEAD iteration (2.1.1)	single, w.r.t. cur weight (1.2.1)
	MOEAD-ilter3	weight-centered clustering [9] (3.3) 1 model per weight (1.2)	1 MOEAD iteration (2.1.1)	single, w.r.t. cur weight (1.2.1)
	MOEAD-ilter4	full dataset (1)	1 model per weight (1.2)	1 MOEAD iteration (2.1.1)	single, w.r.t. cur weight (1.2.1)
	MOEAD-RBF [25]	best per weight (3.2)	3 models (kernels) per objective (2.1) 1 full MOEAD on model response (2.1.2) batch, w.r.t. cur weights (2.2.2)
	MOEAD-EGO [27] fuzzy clustering (3.1)	1 model per cluster and objective (1.1) 1 full MOEAD on EI (2.2.2)	batch, w.r.t. cur weights (2.2.2)
	M-EGO [9]	weight-centered clustering [9] (3.3) 1 model per weight (2.2)	1 full scalar EA on EI (1.2.1)	single, w.r.t. cur weight (1.2.1)
	ParEGO [12]	best sol. for considered weight (3.3) 1 model for 1 weight (1.2)	1 full scalar EA on EI (1.2.1)	single, w.r.t. cur weight (1.2.1)

Table 2 :

 2 12 considered bi-objective benchmark functions. GB of RAM. We rely on the scikit-learn library[START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. while matching the coniguration described in the original paper of each algorithm. For iltering approaches, we use a support vector regression (SVR) with a radial basis function (RBF) kernel. The size of the training set is bounded by 100 when clustering is applicable. It is made of all solutions evaluated so far otherwise.

		objective 1	objective 2
				f 1 f 8 f 14 f 15 f 20
	f 1	Sphere	separable	✓ ✓ ✓	✓	✓
	f 8	Rosenbrock	moderate	✓ ✓	✓	✓
	f 14 Sum of powers ill-conditioned	✓
	f 15 Rastrigin	multi-modal		✓
	f 20 Schwefel	weakly-structured			✓
	a number of function evaluations: ranging from a tight budget of
	500 evaluations to a larger budget of 2 500. For MOEAD-RBF and
	MOEAD-EGO, the number of evaluations per generation is 10. All
	algorithms are implemented in Python, and the experiments are
	conducted on a Intel(R) Core(TM) i5-5200U CPU with 2.20 GHz
	and 8

Table 3 :

 3 Comparison of the competing algorithms with respect to the hypervolume relative deviation under a local reference point for 500 (left) and 2 500 (right) calls to the evaluation function. The rank stands for the number of algorithms that statistically outperform the one under consideration w.r.t a Mann-Whitney test with a p-value of 0.05 and a Bonferroni correction (lower is better). Bold values correspond to the best-performing algorithm for the problem under consideration.

			RS	LHS	MOEAD	MOEAD-ilter1	MOEAD-ilter2	MOEAD-ilter3	MOEAD-ilter4	MOEAD-RBF	MOEAD-EGO	M-EGO	ParEGO	RS	LHS	MOEAD	MOEAD-ilter1	MOEAD-ilter2	MOEAD-ilter3	MOEAD-ilter4	MOEAD-RBF	MOEAD-EGO	M-EGO	ParEGO
							#eval = 500									#eval = 2500
	d = 2	f1-f1 f1-f8 f1-f14 f1-f15 f1-f20 f8-f8 f8-f14 f8-f15 f8-f20 f14-f14 f15-f15 f20-f20	4 6 5 4 6 7 6 7 7 6 6 7	6 4 7 4 4 8 7 7 5 9 5 4	4 2 3 3 3 4 1 4 3 4 3 4	4 3 3 3 0 0 0 0 0 3 3 0	4 1 3 0 3 0 1 0 0 3 1 1	1 2 0 1 0 0 0 0 0 0 1 0	4 0 3 0 1 0 1 0 0 0 0 0	0 0 0 1 0 0 0 0 0 0 2 0	1 0 0 0 0 0 0 0 0 0 0 0	1 1 0 0 1 0 1 0 0 0 0 0	1 4 3 1 0 0 1 0 1 0 0 1	8 7 8 7 8 8 6 9 8 8 9 9	8 7 8 7 7 8 6 7 7 8 9 9	3 0 1 0 0 0 0 0 0 2 2 2	3 2 5 1 0 0 0 0 0 2 2 1	0 0 0 1 0 0 2 0 0 0 2 0	0 0 0 0 0 0 0 0 0 0 0 0	1 0 0 0 0 0 0 0 0 0 0 0	0 0 1 1 0 0 0 0 0 2 2 2	3 0 0 0 0 0 0 0 0 4 0 1	3 5 1 3 5 0 3 0 4 2 2 5	4 7 6 4 5 6 4 0 6 4 7 6
	d = 3	f1-f1 f1-f8 f1-f14 f1-f15 f1-f20 f8-f8 f8-f14 f8-f15 f8-f20 f14-f14 f15-f15 f20-f20	6 6 4 7 6 6 9 6 4 8 4 8	5 5 5 8 9 8 9 8 9 8 4 8	4 1 3 3 4 3 2 5 3 3 4 4	3 0 3 1 1 1 3 3 2 1 4 4	4 1 1 2 3 3 1 3 1 1 4 3	3 1 1 2 3 3 3 0 2 2 4 3	3 1 1 2 2 2 3 3 1 2 4 3	0 0 0 0 0 0 0 0 0 0 0 0	0 0 0 0 0 0 0 0 0 0 0 0	0 1 1 0 1 0 0 0 0 0 0 0	2 1 1 1 3 0 0 0 0 1 0 0	8 9 8 7 8 9 8 9 9 8 9 9	8 5 8 7 8 9 8 9 9 8 8 9	3 1 1 1 1 1 1 1 2 1 2 3	4 0 5 1 2 1 4 0 4 1 4 2	0 2 1 0 1 0 1 0 0 1 0 0	0 0 0 0 0 0 1 0 0 0 0 0	1 2 0 0 0 0 1 0 1 1 0 0	2 0 1 1 0 1 0 1 0 1 0 2	3 0 1 3 0 3 0 2 0 1 0 3	3 2 2 3 5 1 1 1 3 2 0 2	6 3 5 6 7 3 3 4 4 6 5 5
	d = 5	f1-f1 f1-f8 f1-f14 f1-f15 f1-f20 f8-f8 f8-f14 f8-f15 f8-f20 f14-f14 f15-f15 f20-f20	7 9 9 5 9 8 9 8 8 9 8 8	7 8 8 6 9 8 9 8 8 9 8 9	6 4 5 5 5 8 4 8 4 3 3 8	4 4 3 3 3 4 3 3 3 3 3 3	4 4 3 4 3 4 3 3 3 3 3 3	4 4 3 4 3 4 3 3 3 3 3 3	4 4 3 3 3 4 3 3 3 3 3 3	0 0 0 0 0 1 0 0 0 0 0 0	1 0 0 0 0 1 0 0 0 3 0 2	1 0 0 0 1 0 0 0 0 0 0 0	2 2 2 1 2 0 0 0 0 0 0 1	9 9 9 8 9 9 9 9 8 8 8 9	9 9 9 8 9 9 8 9 8 9 8 9	4 1 0 1 1 1 1 4 2 1 2 1	4 1 5 5 1 1 0 4 1 6 3 7	1 3 0 1 1 0 2 0 0 2 0 0	0 1 0 0 1 0 1 0 0 0 0 0	1 1 0 1 1 0 1 0 0 1 0 0	1 0 0 0 0 0 0 0 0 0 0 1	4 1 0 0 1 1 0 0 2 7 0 3	4 1 0 1 2 0 1 1 1 1 2 1	6 5 4 1 5 0 1 3 4 2 2 4
	d = 10	f1-f1 f1-f8 f1-f14 f1-f15 f1-f20 f8-f8 f8-f14 f8-f15 f8-f20 f14-f14 f15-f15 f20-f20	9 9 9 8 9 9 9 9 9 9 6 9	9 9 9 8 9 9 9 9 9 9 7 9	7 7 3 5 6 8 7 7 8 4 3 7	4 4 3 3 4 4 3 3 3 3 3 4	4 4 3 3 4 4 3 3 3 3 3 4	4 4 3 3 4 4 3 3 3 3 3 4	4 4 3 3 4 4 3 3 3 3 3 4	0 0 0 0 0 1 0 0 2 2 2 0	2 1 1 1 1 2 0 0 2 0 1 2	1 1 1 1 1 1 0 0 1 0 0 2	1 1 1 1 1 0 0 0 0 0 0 0	9 9 9 9 9 9 9 9 9 9 9 9	9 9 9 9 9 9 9 9 9 9 9 9	4 3 3 3 4 5 3 3 3 3 3 4	7 3 4 6 8 7 5 3 5 3 6 8	4 3 4 3 1 4 3 3 3 3 3 4	2 3 4 3 4 3 3 3 3 3 3 2	2 3 2 3 4 1 3 3 3 3 3 4	0 0 0 0 0 0 0 0 0 0 0 0	2 1 1 1 1 1 1 1 2 1 2 2	2 1 2 2 2 2 1 1 2 1 1 2	1 1 1 1 1 1 1 1 1 0 0 1
		avg	7.2 7.5 4.4 2.6 2.6 2.4 2.4 0.2 0.4 0.4 0.8 8.5 8.3 1.8 3.1 1.2 0.8	1	0.4 1.2 1.9 3.4

ACKNOWLEDGMENTS

This work was supported by the French national research agency (ANR-16-CE23-0013-01) and the Research Grants Council of Hong Kong (RGC Project No. A-CityU101/16).