Conceptio › Archive › arXiv CS

arXiv CSopen access

Do We Really Need to Approach the Entire Pareto Front in Many-Objective Bayesian Optimisation?

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

artificialintelligence knowledgerepresentation reasoning

artificial intelligence, reasoning, knowledge representation

D O W E R EALLY N EED TO A PPROACH THE E NTIRE PARETO F RONT IN M ANY-O BJECTIVE BAYESIAN O PTIMISATION ?

arXiv:2604.09417v1 [cs.AI] 10 Apr 2026

A P REPRINT Chao Jiang School of Computer Science University of Birmingham [email protected]

Jingyu Huang School of Mathematics University of Birmingham [email protected]

Miqing Li∗ School of Computer Science University of Birmingham [email protected]

April 13, 2026

A BSTRACT Many-objective optimisation, a subset of multi-objective optimisation, involves optimisation problems with more than three objectives. As the number of objectives increases, the number of solutions needed to adequately represent the entire Pareto front typically grows substantially. This makes it challenging, if not infeasible, to design a search algorithm capable of effectively exploring the entire Pareto front. This difficulty is particularly acute in the Bayesian optimisation paradigm, where sample efficiency is critical and only a limited number of solutions (often a few hundred) are evaluated. Moreover, after the optimisation process, the decision-maker eventually selects just one solution for deployment, regardless of how many high-quality, diverse solutions are available. In light of this, we argue an idea that under a limited evaluation budget, it may be more useful to focus on finding a single solution of the highest possible quality for the decision-maker, rather than aiming to approximate the entire Pareto front as existing many-/multi-objective Bayesian optimisation methods typically do. Bearing this idea in mind, this paper proposes a single point-based multi-objective search framework (SPMO) that aims to improve the quality of solutions along a direction that leads to a good tradeoff between objectives. Within SPMO, we present a simple acquisition function, called expected singlepoint improvement (ESPI), working under both noiseless and noisy scenarios. We show that ESPI can be optimised effectively with gradient-based methods via the sample average approximation (SAA) approach and theoretically prove its convergence guarantees under the SAA. We also empirically demonstrate that the proposed SPMO is computationally tractable and outperforms state-of-the-arts on a wide range of benchmark and real-world problems.

Introduction

Multi-objective optimisation problems (MOPs) [Emmerich and Deutz, 2018, Zheng and Wang, 2024] involve scenarios where multiple objectives need to be optimised simultaneously. Unlike single-objective optimisation problems which typically have a single optimal solution, in MOPs there is a set of optimal solutions known as Pareto optimal solutions. The corresponding points in the objective space form what is known as the Pareto front. In general, a multi-objective optimisation algorithm aims to generate a set of solutions that well approximate the Pareto front, from which the decision-maker chooses a solution to deploy based on their preferences. In modern applications, optimisation is becoming increasingly complex, with a growing number of requirements and objectives that need to be considered at the same time [Matrosov et al., 2015, Hierons et al., 2020, Lin et al., 2025]. Taking the car cab design as an example, there are up to nine objectives to be optimised [Deb and Jain, 2013], including cabin space, fuel efficiency, acceleration time, and road noise at various speeds. This has given rise to a new research topic, many-objective optimisation, focusing on MOPs involving more than three objectives [Ishibuchi et al., 2008, Li et al., 2015a]. ∗

Corresponding author: [email protected]

arXiv Template

A P REPRINT

At the same time, many real-world multi-/many-objective optimisation problems are black-box and costly in terms of solution evaluation. This is evident across a range of fields, including chemistry [Park et al., 2018, Shields et al., 2021, Dunlap et al., 2023], materials science [Liang et al., 2021, Low et al., 2024, Peng et al., 2025], and transportation [Deb and Jain, 2013, Jain and Deb, 2013a, Cheaitou and Cariou, 2019, Deb et al., 2009]. For example, in vehicle design optimisation, it can take about 20 hours to evaluate a vehicle design [Youn et al., 2004, Daulton et al., 2021]. To tackle such problems, multi-objective Bayesian optimisation (MOBO) is a very effective approach [Garnett, 2023], along with other alternatives (e.g., surrogate-assisted evolutionary algorithms [Jin, 2011, Liang et al., 2024]). Over the past decades, a variety of effective MOBO methods have emerged, including scalarisation-based methods [Knowles, 2006, Paria et al., 2020, Lin et al., 2022] which convert a multi-objective problem into a number of single-objective problems, and Pareto-based methods which consider Pareto dominance relations over objectives [Daulton et al., 2020, Emmerich et al., 2006, Tu et al., 2022]. Most works aim to find a good approximation of the entire Pareto front. However, with the increase of the objective number, the number of solutions needed to adequately represent the problem’s Pareto front typically grows substantially. For example, for a 10-objective problem, it normally needs 220 points (i.e., 12 3 ) even if only three divisions on each objective are considered [Das and Dennis, 1998]. This difficulty is especially pronounced in Bayesian optimisation, where often only a few hundred solutions can be generated and evaluated. With such a limited budget, it is highly unlikely for an optimisation algorithm to reach or even be close to the Pareto front. On top of that, after the optimisation process, the decision-maker eventually selects just one solution to deploy, regardless of how many high-quality, diverse solutions are provided. Given the above, this paper argues an idea that under a very limited evaluation budget, it may be more useful to focus on finding a single solution of the highest possible quality for the decision-maker, rather than aiming to approximate the entire Pareto front. That is, we may not need to care about diversifying solutions to represent the entire Pareto front, but focus on improving the quality of a single solution. Figure 1 illustrates this idea in a bi-objective case. As can be seen from the figure, in contrast to aiming for a set of diversified solutions which existing MOBO methods typically do [Li et al., 2025, Lin et al., 2022], our method aims for a single solution with better convergence (i.e., closer to the Pareto front). Although it may yield a worse hypervolume (HV) value [Zitzler and Thiele, 1999] compared to the diverse solution set obtained by existing methods, it may be more likely to be chosen by the decision-maker as it achieves a more favourable trade-off among the objectives.

Pareto front Reference point of HV Solution typically obtained by our method Solutions typically obtained by existing methods

Figure 1: An illustration of our idea in a bi-objective case, in comparison with existing methods that aim to search for the entire Pareto front. The red points represent (nondominated) solutions that existing methods may obtain, and the blue point represents what our method aims for. It can be seen that the red points are much more diversified, hence, as a whole, having a better hypervolume (HV) value [Zitzler and Thiele, 1999]. However, the blue point has better convergence (i.e., closer to the Pareto front) than any single red point, which may be more likely to be preferred by the decision-maker.

f1 Bearing this idea in mind, this paper proposes a single point-based multi-objective search framework, called SPMO. The contributions of this work can be summarised as follows. • We propose a novel MOBO search framework that does not aim to approximate the entire Pareto front, but rather focuses on improving the quality of solutions along a single direction that leads to a good tradeoff between objectives. • Within SPMO, we present a simple acquisition function, called Expected Single-Point Improvement (ESPI). We show that ESPI can be optimised effectively with gradient-based methods via the sample average approximation (SAA) approach from Balandat et al. [2020] and also theoretically prove its theoretical convergence guarantees under the SAA. • We consider both noiseless and noisy cases, resulting in two versions of the proposed ESPI. • We verify SPMO through an extensive experimental study, including in comparison with various state-of-thearts, under both sequential and batch optimisation settings, through sensitivity analysis, with different metrics embedded, and on a range of benchmark and real-world problems. 2

arXiv Template

Background and Related Work

2.1

Background

A P REPRINT

Many-Objective Optimisation. Many-objective optimisation, a subset of multi-objective optimisation, refers to an optimisation scenario having more than three objectives to be considered. Without loss of generality, this paper considers the problem of minimising a vector-valued function: f (x) : X → Rm , where x ∈ X (X ⊂ Rd ) and m is the number of objectives. In multi-/many-objective optimisation, a solution x1 is said to dominate x2 , denoted by x1 ≺ x2 , if ∀i ∈ {1, ..., m}, fi (x1 ) ≤ fi (x2 ) and ∃j ∈ {1, ..., m}, fj (x1 ) < fj (x2 ). If a solution x1 ∈ X is not dominated by any other solution, then x1 is said to be Pareto optimal. The collection of Pareto optimal solutions of a problem is called the Pareto set, and its mapping to the objective space is called Pareto front. Bayesian Optimisation (BO). BO is a sample-efficient global optimisation approach that builds a probabilistic surrogate, typically a Gaussian process (GP), and uses an acquisition function α(x) : X → R to decide which points to evaluate. In this work, we model each objective with an independent Gaussian process fi ∼ GP(mi (x), ki (x, x′ )), where mi (x) : X → R is the ith mean function, and ki (·, ·) : X × X → R is the ith covariance function. Given n observed points Dn = {(xt , y t )}nt=1 where y t = f (xt ) + ζ t and the noise ζ t ∼ N (0, diag(σζ2 )), the posterior distribution of the ith objective at a new location x is a Gaussian distribution: p(fi (x)|Dn ) ∼ N (µi (x), σi2 (x)) where µi (x) and σi2 (x) are the mean and variance at x, respectively. Detailed expressions of the mean and variance are given in Appendix A. 2.2

Related Work

Over the past decades, various MOBO methods have been proposed [Konakovic Lukovic et al., 2020, Daulton et al., 2022a]. They can be loosely divided into scalarisation-based and Pareto-based methods. In scalarisation-based methods [Knowles, 2006, Paria et al., 2020], a multi-objective problem is converted into a number of single-objective problems [Chugh, 2020]. Hence, one can leverage acquisition functions [Lai and Robbins, 1985] from single-objective BO to decide which point to evaluate. For instance, using random augmented Tchebycheff scalarisations [Miettinen, 1999], ParEGO [Knowles, 2006] and TS-TCH [Paria et al., 2020] optimise expected improvement (EI) [Jones et al., 1998] and Thompson sampling (TS) [Thompson, 1933], respectively. In contrast, Pareto-based methods consider Pareto dominance relations over objectives [Emmerich et al., 2006, Tu et al., 2022]. A popular idea is to use HV as maximising the HV value is equivalent to finding the entire Pareto front [Shang et al., 2020]. Expected hypervolume improvement (EHVI) is commonly used in MOBO [Couckuyt et al., 2014, Daulton et al., 2020, 2021]. Another idea is to leverage information theory to guide exploration toward regions likely contributing to the Pareto front. For instance, joint entropy search (JES) [Tu et al., 2022] selects points that maximise the joint information gain for optimal inputs (i.e., the approximated Pareto set) and outputs (i.e., the approximated Pareto front). All of the above methods aim to approach the entire Pareto front. It is worth noting that, similar to our approach, a few studies do not attempt to approximate the entire Pareto front. Some methods instead target a specific region of the front, such as the central area [Gaudrie et al., 2018, 2020, Binois et al., 2020]. Another line of work incorporates decision-maker preferences by dynamically adjusting the target region based on elicited or updated preferences during optimisation [Abdolshah et al., 2019, Astudillo and Frazier, 2020, Ozaki et al., 2024, Ip et al., 2025]. In contrast, our method assumes no prior knowledge of decision-maker preferences and seeks to identify a high-quality trade-off solution across objectives. A more detailed discussion of related work can be found in Appendix C.

The Proposed Method

In this section, we first give the proposed MOBO framework. We then present the considered acquisition function (called ESPI), which is based on a simple distance-based metric. We note that analytically solving ESPI is not feasible, and thus consider its Monte Carlo approximation. Lastly, we consider ESPI under noisy cases, namely noisy ESPI (NESPI), and also its MC approximation. Single Point-based Multi-Objective (SPMO) Framework. Algorithm 1 gives the procedure of the proposed SPMO framework. As can be seen, SPMO is very similar to a standard MOBO algorithm, except for the step of maximising the acquisition function based on a single-point quality metric (line 3). In principle, any metric that can reflect the quality of a solution in achieving a good trade-off between objectives can be adopted. This includes distance-based metrics and scalarisation-based metrics, such as the weighted sum or augmented Tchebycheff [Miettinen, 1999] with 3

arXiv Template

A P REPRINT

Algorithm 1: Single Point-based Multi-Objective (SPMO) Framework Input: f : Expensive black-box problem with m objectives; T : Maximum number of evaluations; g: Metric that measures the quality of a single point; α: Acquisition function; 0 : Initial observed points. Dn0 := {(xt , y t )}nt=1 1 for n = n0 + 1 : T do 2 GP s ← Train GPs(Dn ) // Train m Gaussian process models 3 xn ← arg maxx∈X α g(x, GP s) // Maximise the acquisition function based on the single-point quality metric g(·) 4 y n ← f (xn ) + ζ n // Evaluate the solution xn n n−1 n n 5 D ←D ∪ {(x , y )} // Augment the observed solution Output: DT : Observed solutions.

1 1 a fixed weight vector (e.g., ( m ,..., m ) ∈ Rm in the m-objective case). Here, we consider a simple distance-based metric, and we will compare different metrics in our experiments (Section 6; details in Appendix F.5).

Single-Point We consider the distance of solutions to a utopian point: g(f (x), z ∗ ) = ∥f (x) − q Improvement (SPI). P m ∗ 2 , where m denotes the number of objectives, and z ∗ = (z ∗ , z ∗ , . . . , z ∗ ) is a utopian z∗∥ = m 1 2 i=1 fi (x) − zi point, i.e., zi∗ ≤ minx∈X fi (x). In many real-world cases, the utopian value of an objective can be loosely estimated, for example, by assuming idealised conditions such as zero cost, time or error [Branke et al., 2008]. In our experimental evaluation, we perform a sensitivity analysis on the choice of the utopian point, and it shows that substantially different settings can yield consistent results. We first consider noiseless cases, i.e., ȳ t = f (x̄t ). Let D̄n = {(x̄t , ȳ t }nt=1 be n observed points and X̄ n = {x̄t }nt=1 be the set of all the observed decision vectors in D̄n . For any point x ∈ X , we define the single-point improvement (SPI) as: ISP (f (x)|g ∗ , z ∗ , D̄n ) = max 0, g ∗ − ∥f (x) − z ∗ ∥ (1) where g ∗ = minx∈X̄ n g(f (x), z ∗ ). Expected Single-Point Improvement (ESPI). We now present ESPI to account for the posterior distribution p(f |D̄n ). Suppose that we independently model each objective fi as a Gaussian process based on D̄n , Then, the posterior of each fi at a new location x is a Gaussian random variable, i.e., p(fi (x)|D̄n ) ∼ N µi (x), σi2 (x) , in which f1 , . . . , fm are mutually independent Gaussians. Let ηi := fi − zi∗ . Then we obtain ηi ∼ N µi (x) − zi∗ , σi2 (x) which is a Gaussian as well. The proposed ESPI is defined as: αESPI (x) = E ISP (f (x)|g ∗ , z ∗ , D̄n ) (2) = Ep(η) max(0, g ∗ − ∥η∥)] An illustration of ESPI in a bi-objective case is given in Figure 5 of Appendix B for aiding understanding. Note that the integral in Eq. 2 cannot be solved analytically as it involves the distribution of ∥η∥, whose PDF and CDF have no closed-form expressions and are typically computed via numerical methods [Imhof, 1961, Ruben, 1962, Das, 2025]. Hence, we use the MC integration with samples from the posterior f˜t (x) ∼ p(f (x)|D̄n ) for t = 1, . . . , N to estimate Eq. 2: N 1 X αESPI (x) ≈ α̂ESPI (x) = ISP (f˜t (x)|g ∗ , z ∗ , D̄n ) (3) N t=1 Noisy Expected Single-Point Improvement (NESPI). In the real world, it is not uncommon to encounter an optimisation problem with noises: y t = f (xt ) + ζ t , where ζ t ∼ N (0, diag(σζ2 )). In noisy cases, simply using the observed best distance g ∗ may adversely affect the optimisation performance. Here, we present an extension of ESPI, i.e., noisy ESPI (NESPI). Let Dn = {(xt , y t )}nt=1 be n observed points and X n = {xt }nt=1 be the set of all the observed decision vectors in Dn . By considering the uncertainty in the function values at X n , the proposed NESPI is defined as: Z αNESPI (x) = αESPI (x|gˆ∗ )p(f |Dn )df (4) 4

arXiv Template

A P REPRINT

where gˆ∗ denotes the smallest distance to the utopian point over f (X n ). Note that in noiseless cases, NESPI is equivalent to ESPI. Additionally, ESPI and NESPI can be naturally extended to the parallel (batch) setting by using sequential greedy approximation [Balandat et al., 2020]. Like in the noiseless case, the integral in Eq. 4 is also analytically intractable but can be approximated using the MC integration. Let f˜t (x) ∼ p(f (x)|Dn ) for t = 1, . . . , N be samples from the posterior, and let gˆ∗ = minx∈X n ∥f˜t (x)− z ∗ ∥ be the smallest distance to utopian point over the previously evaluated points under the sampled function f˜t (x). PN Then, αNESPI ≈ N1 t=1 αESPI (x|gˆ∗ , z ∗ , Dn ). Using the MC integration, the inner expectation in αNESPI can be computed simultaneously using samples from the joint posterior f˜t (x, X n ) ∼ p(f (x, X n )|Dn ) over x and X n : N

αNESPI (x) ≈ α̂NESPI (x) =

1 X ISP (f˜t |gˆ∗ , z ∗ , Dn ) N t=1

(5)

Optimising ESPI and NESPI

Having presented the MC estimators of ESPI and NESPI, we are now ready to optimise them. Differentiability. The MC estimators of ESPI and NESPI (α̂ESPI (x) in Eq. 3 and α̂NESPI (x) in Eq. 5) are differentiable with respect to x. We are able to automatically compute exact gradients of the MC estimators of ESPI and NESPI (∇x α̂ESPI (x) and ∇x α̂NESPI (x)) by leveraging the auto-differentiation in modern computational frameworks. This facilitates efficient gradient-based optimisation of ESPI and NESPI. SAA Convergence Results. The sample average approximation (SAA) approach [Kleywegt et al., 2002], which addresses stochastic optimisation problems by using the MC simulation, has gained increasing popularity and has become a standard technique in BO for optimising MC-based acquisition functions [Balandat et al., 2020]. By fixing the base samples, the SAA yields a deterministic acquisition function which enables using (quasi-) higher-order optimisation algorithms to obtain fast convergence rates for acquisition optimisation. We now give the theoretical convergence guarantees of ESPI under the SAA. Theorem 4.1. Suppose that X is compact and f has a multi-output GP prior whose mean and covariance functions are continuously differentiable. Let α∗ESPI := maxx∈X αESPI (x) denote the maximum of ESPI, S ∗ := arg maxx∈X αESPI (x) denote the set of maximisers of αESPI , α̂N ESPI (x) denote the deterministic function via the ∗ N base samples {ϵt }N ∼ N (0, I ). Suppose x̂ ∈ arg max α̂ m x∈X ESPI (x), then t=1 N ∗ ∗ (1) α̂N ESPI (x̂N ) → αESPI a.s.

(2) d(x̂N∗ , S ∗ ) → 0 a.s., where d(x̂N∗ , S ∗ ) := inf x∈S ∗ ∥x̂N∗ − x∥. The proof of the theorem is given in Appendix D.1. This theorem indicates that one is able to optimise the MC estimator of the acquisition function ESPI to obtain a solution that converges almost surely to the optimal solution of the original function. For the noisy case NESPI, the theorem of the theoretical convergence guarantees under the SAA (Theorem D.1), together with its proof, is provided in Appendix D.2.

Experimental Design

Compared Methods. To evaluate the proposed SPMO, we consider six MOBO methods. They include one baseline method (Sobol [Sobol, 1967]), four well-established methods that aim to approximate the entire Pareto front, and one method that aims at the trade-off region of the Pareto front. The four methods consist of two scalarisation-based methods, ParEGO [Knowles, 2006] (along with its noisy variant NParEGO [Daulton et al., 2021]) and TS-TCH [Paria et al., 2020], and two Pareto-based methods, EHVI [Daulton et al., 2020] (along with its noisy variant NEHVI [Daulton et al., 2021]), and JES [Tu et al., 2022]. For the method that does not aim at the entire Pareto front, we consider C-EHVI [Gaudrie et al., 2018, 2020]. C-EHVI prefers the central region of the Pareto front, and we would like to see if it is competitive against our method in identifying a well-balanced solution. For all EI-based methods, i.e., ParEGO, NParEGO, EHVI, NEHVI, C-EHVI and SPMO, we use the log version as suggested by Ament et al. [2023]. Note that we only consider NESPI in this work, as it is equivalent to ESPI under noiseless cases.2 2

Although NEHVI is equivalent to EHVI under noiseless cases, the wall time of NEHVI is much higher than EHVI (see Table 23). Hence we consider EHVI and NEHVI under noiseless and noisy cases, respectively.

arXiv Template

A P REPRINT

Benchmarks and Real-World Problems. For benchmark problems, we first choose two most widely scalable functions, DTLZ1 and DTLZ2 [Deb et al., 2005]. However, their Pareto fronts are rather homogeneous, i.e., with a simplex shape. We then include their inverted versions, i.e., inverted DTLZ1 and inverted DTLZ2 [Deb and Jain, 2013]. These problems do not include the one with convex Pareto fronts nor different objective scales. We thus add convex DTLZ2 and scaled DTLZ2 [Jain and Deb, 2013a]. We also give the results of other DTLZ problems which have different features (e.g., degenerate and disconnected Pareto fronts) [Deb et al., 2005, Cheng et al., 2017] in Appendix F.1; those problems are widely used in many-objective optimisation [Deb and Jain, 2013, Li et al., 2014a,b, 2015b]. Each problem is considered with 3, 5 and 10 objectives, following the practice in Deb and Jain [2013], Jain and Deb [2013a], Li et al. [2014b], and tested under both noiseless and noisy cases. We also consider two well-studied expensive real-world problems [Tanabe and Ishibuchi, 2020], i.e., car side impact design [Jain and Deb, 2013a] and car cab design [Deb and Jain, 2013]. The former is a four-objective problem without noise. The latter, which has nine objectives, involves four stochastic variables (out of the total seven variables) that introduce noise into the optimisation process, thus a natural optimisation problem with noise. For the other problems without noise, to make their noise cases, we use the additive zero-mean Gaussian noise with a standard deviation of 0.1, as suggested in Hernandez-Lobato et al. [2016], Jiang and Li [2025a]. The details of the problem formulations are given in Appendix E.3. Performance Metrics. The proposed method aims to find a single trade-off solution between objectives which has a high chance of being favoured by the decision-maker, and therefore requires the use of appropriate metrics to fairly evaluate our method [Li and Yao, 2019, Li et al., 2022]. As such, we first consider two single-point-based metrics, i.e., the distance-based metric used in the proposed method (reported as log-distance for better visualisation) and the HV-based metric [Zitzler and Thiele, 1999], which measures the HV contribution of a solution. We expect that our method performs well on the distance-based metric as we directly optimise it. Notably, we are not certain whether our method performs best on the single-point HV since we do not directly optimise it. On the other hand, since most of the compared methods aim to achieve a good approximation of the entire Pareto front, we also consider the HV of the whole (nondominated) solution set obtained. It is expected that our method performs poorly, compared to other methods since we only optimise one point, rather than maximising the HV of the whole set (see Figure 1). For the reference point of the two HV metrics, we followed the practice in Ishibuchi et al. [2018], Balandat et al. [2020], Chugh [2020], Daulton et al. [2020] (see detailed settings in Appendix E.2). Budget and Statistical Validation. For all the methods, we allow a maximum of 200 evaluations, following the practice in Daulton et al. [2020, 2021], Konakovic Lukovic et al. [2020]. To enable statistical comparisons, each optimisation was repeated 30 times. We use the Wilcoxon rank-sum test [Wilcoxon, 1992] at a significance level of α = 0.05 and Holm-Bonferroni correction [Holm, 1979] to see if our method differs significantly from each peer method.

Experimental Results

We first report the results under noiseless cases, then under noisy cases and under batch settings. Next, we perform the sensitivity analysis of the utopian point used. Then, we compare the proposed framework working with different singlepoint metrics, e.g., the distance-based, weighted sum-based and Tchebycheff-based. Lastly, we give the acquisition optimisation wall time for all the algorithms. Noiseless Cases. We begin our evaluation by considering the distance metric, for which we expect good performance obtained by our SPMO. Table 1 shows the results of SPMO and the other methods on the seven benchmark (with five objectives) and real-world problems. Unsurprisingly, as can be seen from the table, our method significantly outperforms the other methods on all the problems. In addition, to understand the anytime performance, Figure 2 presents the trajectories of the distance-based metric obtained by the six methods. As seen, SPMO demonstrates a clear advantage over the other methods, obtaining a better convergence rate from the very beginning. We now compare the methods by the HV metric of a single solution. Table 2 shows the results of the best solution (in terms of its HV value) obtained by SPMO and the five peer methods. As can be seen, SPMO is still very competitive, performing best on all the problems except DTLZ2 where EHVI obtains the best HV. To get a sense of what such a best-HV solution looks like, we use a spider chart to plot the solution on the five-objective inverted DTLZ1 problem in Figure 3. As seen, the solution of SPMO has the largest area, and it actually Pareto dominates the solutions of the other methods (i.e., better or at least equal on all five objectives). Lastly, we consider the HV results of all the solutions obtained by the compared methods (Table 3). Interestingly, although SPMO does not aim to approximate the entire Pareto front, it still gets fairly good results, outperforming the other methods on at least 3 out of the 7 problems. One explanation for this is that within very tight budget, searching 6

arXiv Template

A P REPRINT

Table 1: Results of the distance-based metric (log distance) obtained by the SPMO and the six peer methods on the benchmark problems with 5 objectives and the car side impact design problem on 30 runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Car side impact Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI C-EHVI JES SPMO

3.7e+0 (3.0e–1)+ 3.4e+0 (1.9e–1)+ 3.9e+0 (2.2e–1)+ 3.5e+0 (9.6e–2)+ 3.6e+0 (1.4e–1)+ 3.4e+0 (1.3e–1)+ 3.1e+0 (3.0e–1)

2.4e–1 (4.6e–2)+ 6.7e–2 (8.1e–2)+ 2.0e–1 (4.1e–2)+ 9.3e–3 (3.5e–3)+ 3.3e–3 (3.8e–3)+ 1.1e–1 (8.0e–2)+ 9.0e–4 (8.3e–4)

4.8e+0 (3.1e–1)+ 3.2e+0 (5.4e–1)+ 4.8e+0 (3.2e–1)+ 4.0e+0 (4.4e–1)+ 3.9e+0 (4.6e–1)+ 4.5e+0 (1.7e–1)+ 2.9e+0 (4.9e–1)

6.0e–1 (5.7e–2)+ 2.5e–1 (1.3e–2)+ 4.6e–1 (1.2e–2)+ 2.3e–1 (5.5e–3)+ 2.5e–1 (1.8e–2)+ 2.6e–1 (2.0e–2)+ 2.1e–1 (5.0e–5)

-3.1e–1 (2.5e–1)+ -1.7e+0 (2.0e–1)+ -7.1e–1 (2.1e–1)+ -1.4e+0 (2.1e–1)+ -3.5e–1 (2.9e–1)+ -1.0e+0 (4.5e–1)+ -2.1e+0 (7.7e–3)

2.3e–1 (5.1e–2)+ 1.3e–1 (8.9e–2)+ 2.6e–1 (3.8e–2)+ 3.1e–1 (6.0e–2)+ 7.6e–2 (8.0e–2)+ 8.8e–2 (9.3e–2)+ 1.7e–4 (1.2e–4)

-1.9e–1 (2.5e–2)+ -3.3e–1 (6.0e–3)+ -2.7e–1 (1.5e–2)+ -3.3e–1 (1.0e–2)+ -3.3e–1 (1.0e–2)+ -3.3e–1 (8.3e–3)+ -3.4e–1 (2.3e–6)

7/ 0/ 0 7/ 0/ 0 7/ 0/ 0 7/ 0/ 0 7/ 0/ 0 7/ 0/ 0

DTLZ1 (5 obj)

4.4

DTLZ2 (5 obj)

DTLZ1 (5)

4.2

Log distance

4.0 3.8 3.6 3.2 25

Log distance

0.6

0.20

4.5

4.2

0.15

4.0

0.05

0.10

100 125 150 175 200

3.8

0.00

3.5 3.0 25

100

0.3 50

100

125

150

175

Number of evaluations

0.25

1.0

0.15

25 25

Sobol

200

100 125 150 175 200

0.20 0.25

0.10

0.30

0.05

5075 10075 125 100 125 15025 175 200 0.35 150 175 200 50 75 100 125 150 175 200 0.00

Number of evaluations Number of evaluations

Number of evaluations

ParEGO

TS-TCH

EHVI

Car side impact design (4 obj)

0.10 0.15

0.20

Number of evaluations

Scaled DTLZ2 (5 obj)

3.4 0.5

200

175

0.30

2.0 25

150

Number of evaluations

3.2 1.5

0.4

125

Convex DTLZ2 (5 obj)

3.6 0.0

0.5

0.2

0.25

Number of evaluations

Inverted DTLZ2 (5 obj) 0.7

5.0

4.4

Log distance

3.4

Inverted DTLZ1 (5 obj)

5.5

0.30

C-EHVI

JES

100 125 150 175 200

Number of evaluations

SPMO (ours)

Figure 2: Trajectories of the distance-based metric (log distance) obtained by the seven methods on the benchmark problems with 5 objectives and the car side impact design problem. Each coloured line represents the mean metric value on 30 independent runs (after the initial Sobol samples, represented by the dashed grey line).

for improving convergence of solutions may play a bigger part than searching for improving diversity, thus contributing more to the HV value. For simplicity, we here only show results of the benchmark problems with five objectives. Results on 3- and 10-objective problems can be found in Appendix F.1. A general pattern is that as the number of objectives increases, the advantages of SPMO become more pronounced. On the 3-objective problems, the differences between SPMO and the peer methods in terms of the two single-point metrics are relatively small, whereas on the 10-objective problems, the gaps in both metrics become substantially larger (see Figures 6 and 7 in the Appendix). Regarding the HV of all evaluated solutions, SPMO statistically outperforms the peer methods on at least 2 and 5 out of the 6 problems on the 3-objective and 10-objective cases, respectively. Noisy Cases. We compare the proposed SPMO with the peer methods on the seven noisy benchmark and realworld problems. The results (the distance-based metric, two HV metrics, and convergence trajectories) are given in Appendix F.2. Like in the noiseless setting, SPMO significantly outperforms the peer methods on all the problems with respect to the single-point metrics (distance and HV). Regarding the HV of all the (nondominated) solutions obtained, SPMO achieves the best performance on the majority of the problems. Figure 4 shows the spider chart of the best solution (with respect to its HV) obtained by each algorithm in a typical run on the 9-objective car cab design problem. The violin plot of their HV values in 30 independent runs is given in the left panel for reference. As can be seen, the solution of SPMO is the best or close to the best on most of the objectives (except the first objective), thus having a clearly larger area. Batch Setting. Previously, we considered the case in the sequential setting, where solutions are evaluated sequentially. Now we want to see if the proposed method works in the batch setting. Here, the batch size q is set to 5, a commonly 7

arXiv Template

A P REPRINT

Table 2: The HV of the best solution (in terms of its HV value) obtained by SPMO and the peer methods on the benchmark problems with 5 objectives and the car side impact design problem on 30 runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Car side impact Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI C-EHVI JES SPMO

8.7e+12 (3.7e+11)+ 9.1e+12 (1.8e+11)+ 8.4e+12 (3.8e+11)+ 9.0e+12 (9.5e+10)+ 9.0e+12 (1.0e+11)+ 9.0e+12 (1.2e+11)+ 9.3e+12 (2.8e+11)

3.5e–2 (1.5e–2)+ 1.3e–1 (5.6e–2)+ 3.2e–2 (1.3e–2)+ 1.9e–1 (7.2e–3)∼ 1.6e–1 (1.9e–2)+ 1.0e–1 (3.9e–2)+ 1.9e–1 (1.4e–2)

5.1e+12 (1.1e+12)+ 8.8e+12 (7.0e+11)+ 4.9e+12 (1.2e+12)+ 7.5e+12 (9.7e+11)+ 7.7e+12 (7.8e+11)+ 6.2e+12 (5.0e+11)+ 9.2e+12 (4.8e+11)

8.7e–4 (1.5e–3)+ 4.0e–2 (2.9e–3)+ 7.1e–3 (1.1e–3)+ 4.5e–2 (1.4e–3)+ 4.1e–2 (4.1e–3)+ 3.9e–2 (4.5e–3)+ 4.9e–2 (1.2e–5)

3.8e–1 (1.7e–1)+ 1.1e+0 (6.2e–2)+ 6.5e–1 (1.4e–1)+ 1.0e+0 (9.7e–2)+ 4.3e–1 (2.1e–1)+ 8.5e–1 (2.4e–1)+ 1.3e+0 (5.0e–3)

3.8e–2 (2.1e–2)+ 8.7e–2 (5.3e–2)+ 2.8e–2 (1.5e–2)+ 1.5e–2 (1.4e–2)+ 1.2e–1 (5.1e–2)+ 1.1e–1 (5.5e–2)+ 1.6e–1 (1.6e–2)

2.6e–1 (1.8e–2)+ 3.6e–1 (2.6e–3)+ 3.1e–1 (9.3e–3)+ 3.6e–1 (6.8e–3)+ 3.6e–1 (4.4e–3)∼ 3.6e–1 (5.1e–3)+ 3.6e–1 (3.2e–3)

7/ 0/ 0 7/ 0/ 0 7/ 0/ 0 6/ 1/ 0 6/ 1/ 0 7/ 0/ 0

Table 3: The HV of all the solutions obtained by the seven methods on the seven benchmark (with five objectives) and real-world problems on 30 independent runs. The method with the best mean HV is highlighted in bold. The symbols “+”, “∼”, and “−” indicate that a method is statistically worse than, equivalent to, and better than SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Car side impact Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI C-EHVI JES SPMO

1.0e+13 (3.1e+10)+ 1.0e+13 (3.8e+10)− 1.0e+13 (2.6e+10)+ 1.0e+13 (3.3e+8)− 1.0e+13 (1.6e+11)∼ 1.0e+13 (9.0e+9)− 1.0e+13 (8.1e+10)

7.9e–2 (2.5e–2)+ 3.6e–1 (2.1e–1)+ 5.9e–2 (2.7e–2)+ 8.3e–1 (6.2e–2)− 3.6e–1 (8.5e–2)+ 2.8e–1 (1.5e–1)+ 5.0e–1 (8.2e–2)

5.6e+12 (8.7e+11)+ 9.1e+12 (6.2e+11)+ 5.4e+12 (1.0e+12)+ 8.2e+12 (6.6e+11)+ 8.0e+12 (6.9e+11)+ 8.6e+12 (3.5e+11)+ 9.6e+12 (2.9e+11)

9.3e–4 (1.5e–3)+ 1.4e–1 (1.0e–2)− 2.0e–2 (3.6e–3)+ 2.1e–1 (1.9e–3)− 7.2e–2 (9.2e–3)− 1.4e–1 (9.4e–3)− 6.4e–2 (2.8e–3)

5.8e–1 (2.1e–1)+ 1.5e+0 (2.8e–2)− 1.0e+0 (1.5e–1)+ 1.5e+0 (6.9e–2)+ 5.8e–1 (2.7e–1)+ 1.3e+0 (2.6e–1)+ 1.5e+0 (3.2e–2)

8.6e–2 (3.1e–2)+ 2.2e–1 (1.6e–1)+ 5.1e–2 (2.7e–2)+ 2.0e–2 (1.7e–2)+ 2.1e–1 (1.1e–1)+ 3.1e–1 (2.0e–1)∼ 3.0e–1 (1.1e–1)

5.2e–1 (1.0e–2)+ 7.4e–1 (1.4e–2)− 6.5e–1 (1.1e–2)− 7.4e–1 (9.6e–3)− 5.9e–1 ( 3.0e–2)− 7.4e–1 (1.3e–2)− 5.3e–1 (3.1e–2)

7/ 0/ 0 3/ 0/ 4 6/ 0/ 1 3/ 0/ 4 4/ 1/ 2 3/ 1/ 3

2 3 100

Sobol ParEGO TS-TCH EHVI C-EHVI JES SPMO

Figure 3: Spider chart of the best solution (in terms of its HV) obtained by the seven methods on the inverted DTLZ1 problem with 5 objectives in a typical run. Each axis in the spider chart represents one objective. Here, the objective values are multiplied by −1 in this minimisation problem, such that a solution with a larger area indicates better quality.

5 used value [Lin et al., 2022]. As can be seen in Tables 14–16 and Figures 13–15 (Appendix F.3), similar to the results in the sequential setting, SPMO generally performs best. On the single-point metrics, it achieves the smallest distance on all the problems and the highest HV on 5 out of the 6 problems (except on DTLZ2). As for the HV of all evaluated solutions, SPMO obtains the best HV on two problems and takes the second or third places on the remaining ones. This indicates that prioritising convergence is also very useful for many-objective problems in batch settings. Sensitivity Analysis. A parameter needed in the proposed method is the utopian point. In our experiment, we set it to be the vector consisting of the best value on each objective (i.e., the problem’s ideal point). However, in real life, the ideal point is usually unknown before the optimisation. Hence, we would like to investigate how much different utopian points affect the performance. In this context, we consider three different settings. The first one is slightly better than the ideal point, i.e., with a difference of 0.01, the second is fairly better than the ideal point (i.e. 0.1), and the last one is significantly better than the ideal point (i.e. 1.0). The results are given in Appendix F.4. As can be seen, interestingly, SPMO with the three lower utopian points with different levels performs better than or at least equivalently to SPMO with our current setting. This indicates that 1) using the ideal point in the proposed method may not be the best choice (though it performs better than the other MOBO methods), and 2) SPMO’s performance is robust to the choice of the utopian point and a liberal estimate is sufficient to achieve good results, e.g., considering zero cost in real-world cases. Comparison of Single-Point Metrics within SPMO. In the proposed SPMO framework, we employ a distance metric (i.e., the distance of a solution to the utopian point), denoted as SPMOdist . However, different metrics can be adopted provided that they can reflect the quality of a solution in achieving a good trade-off between objectives. We now consider two other well-known metrics, weighted sum and Tchebycheff scalarisation (with the same weights 1 1 (m ,..., m )), denoted by SPMOws and SPMOT ch , respectively. We compare these three versions of SPMO. The 8

Car cab design

A P REPRINT

3 2

4.0

4.5 5.0

1.0

5.5

0.8

0.6

0.4

0.0 0.2

SPMO

JES

C-EHVI

TS-TCH

6 NParEGO

6.0

Sobol

Log HV (single point)

arXiv Template

Sobol NParEGO TS-TCH C-EHVI JES SPMO

9 8

Figure 4: Left: Violin plots of the HV values of the best solution (with respect to its HV) obtained by all the methods on the car cab design problem in 30 independent runs. Right: The objective values (normalised and multiplied by −1) of the best solution (with the highest HV) obtained by each method on the cab design problem in a typical run.

results (given in Appendix F.5) show that SPMOdist performs in general better than SPMOws and SPMOT ch . It obtains the best result on at least 4 out of the 6 problems on the HV of the best solution. As for the HV of all the solutions, SPMOdist performs best on DTLZ1 and its variants, but worse than SPMOT ch on DTLZ2 and its variants (except convex DTLZ2). A possible explanation is that SPMOT ch has slower convergence and can be better in exploring different solutions, thus better on relatively easy-to-converge problems. Acquisition Optimisation Wall Time. Lastly, we present the wall time for optimising the acquisition function (i.e., determining a solution to be evaluated). The results (Appendix F.6) show that our method is among the fastest algorithms. When the number of objectives is 3 or 5, the time of all the methods is acceptable with a maximum of 98 seconds. As the number of objectives increases to 10, hypervolume-based methods (i.e., EHVI and NEHVI) become very expensive (taking about half an hour and more than 3 hours, respectively).3 The proposed SPMO method shows high computational efficiency, achieving the lowest time requirement in four out of the six instances.

Conclusion

This work presented a multi-objective BO framework that aims to find a single trade-off solution of the highest possible quality with respect to multiple objectives, rather than seeking to explore their entire Pareto front. We theoretically proved the convergence guarantees under the SAA and empirically verified the proposed framework through extensive experiments, including on noiseless/noisy and sequential/batch cases, by sensitivity analysis, with different metrics for the acquisition function, and on a range of benchmark and real-world problems. A noticeable limitation of the proposed framework is that it focuses on finding a single trade-off point, thus failing to capture the information about the entire Pareto front; it thus may be less useful for certain applications where such information is valuable (e.g., the Pareto front’s ranges and nadir points). A detailed discussion of its applicability is provided in Appendix G. However, notably, the proposed framework showed its competitiveness against existing state-of-the-arts with respect to even the quality of the whole solution set (through HV of all solutions, see Table 3). Future work includes studying and enhancing the scalability of the proposed methods (i.e., in higher-dimensional search space) and extending their applicability to other scenarios, e.g., multi-fidelity optimisation (see Appendix H for more details).

References Michael TM Emmerich and André H Deutz. A tutorial on multiobjective optimization: fundamentals and evolutionary methods. Natural Computing, 17:585–609, 2018. Ruihao Zheng and Zhenkun Wang. Boundary decomposition for nadir objective vector estimation. In Advances in Neural Information Processing Systems, volume 37, pages 14349–14385. Curran Associates, Inc., 2024. Evgenii S Matrosov, Ivana Huskova, Joseph R Kasprzyk, Julien J Harou, Chris Lambert, and Patrick M Reed. Manyobjective optimization and visual analytics reveal key trade-offs for london’s water supply. Journal of Hydrology, 531:1040–1053, 2015. 3

Wall time is measured based on the initial samples. Due to the exponentially increasing computational complexity with objectives, EHVI and NEHVI are fully evaluated only for problems with objectives m ≤ 5.

arXiv Template

A P REPRINT

Robert M Hierons, Miqing Li, Xiaohui Liu, Jose Antonio Parejo, Sergio Segura, and Xin Yao. Many-objective test suite generation for software product lines. ACM Transactions on Software Engineering and Methodology (TOSEM), 29(1):1–46, 2020. Xi Lin, Yilu Liu, Xiaoyuan Zhang, Fei Liu, Zhenkun Wang, and Qingfu Zhang. Few for many: Tchebycheff set scalarization for many-objective optimization. In The Thirteenth International Conference on Learning Representations, 2025. Kalyanmoy Deb and Himanshu Jain. An evolutionary many-objective optimization algorithm using reference-pointbased nondominated sorting approach, part I: solving problems with box constraints. IEEE transactions on Evolutionary Computation, 18(4):577–601, 2013. Hisao Ishibuchi, Noritaka Tsukamoto, and Yusuke Nojima. Evolutionary many-objective optimization: A short review. In 2008 IEEE congress on evolutionary computation (IEEE world congress on computational intelligence), pages 2419–2426. IEEE, 2008. Bingdong Li, Jinlong Li, Ke Tang, and Xin Yao. Many-objective evolutionary algorithms: A survey. ACM Computing Surveys (CSUR), 48(1):1–35, 2015a. Seongeon Park, Jonggeol Na, Minjun Kim, and Jong Min Lee. Multi-objective Bayesian optimization of chemical reactor design using computational fluid dynamics. Computers & Chemical Engineering, 119:25–37, 2018. Benjamin J Shields, Jason Stevens, Jun Li, Marvin Parasram, Farhan Damani, Jesus I Martinez Alvarado, Jacob M Janey, Ryan P Adams, and Abigail G Doyle. Bayesian reaction optimization as a tool for chemical synthesis. Nature, 590(7844):89–96, 2021. John H Dunlap, Jeffrey G Ethier, Amelia A Putnam-Neeb, Sanjay Iyer, Shao-Xiong Lennon Luo, Haosheng Feng, Jose Antonio Garrido Torres, Abigail G Doyle, Timothy M Swager, Richard A Vaia, et al. Continuous flow synthesis of pyridinium salts accelerated by multi-objective Bayesian optimization with active learning. Chemical Science, 14 (30):8061–8069, 2023. Qiaohao Liang, Aldair E Gongora, Zekun Ren, Armi Tiihonen, Zhe Liu, Shijing Sun, James R Deneault, Daniil Bash, Flore Mekki-Berrada, Saif A Khan, et al. Benchmarking the performance of Bayesian optimization across multiple experimental materials science domains. npj Computational Materials, 7(1):188, 2021. Andre KY Low, Flore Mekki-Berrada, Abhishek Gupta, Aleksandr Ostudin, Jiaxun Xie, Eleonore Vissol-Gaudin, YeeFun Lim, Qianxiao Li, Yew Soon Ong, Saif A Khan, et al. Evolution-guided Bayesian optimization for constrained multi-objective optimization in self-driving labs. npj Computational Materials, 10(1):104, 2024. Peng Peng, Yi Peng, Fuguo Liu, Shuai Long, Cheng Zhang, Aitao Tang, Jia She, Jianyue Zhang, and Fusheng Pan. Bayesian optimization and explainable machine learning for high-dimensional multi-objective optimization of biodegradable magnesium alloys. Journal of Materials Science & Technology, 2025. Himanshu Jain and Kalyanmoy Deb. An evolutionary many-objective optimization algorithm using reference-point based nondominated sorting approach, part II: Handling constraints and extending to an adaptive approach. IEEE Transactions on Evolutionary Computation, 18(4):602–622, 2013a. Ali Cheaitou and Pierre Cariou. Greening of maritime transportation: a multi-objective optimization approach. Annals of Operations Research, 273(1):501–525, 2019. Kalyanmoy Deb, Shubham Gupta, David Daum, Jürgen Branke, Abhishek Kumar Mall, and Dhanesh Padmanabhan. Reliability-based optimization using evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 13 (5):1054–1074, 2009. Byeng D Youn, KK Choi, R-J Yang, and Lei Gu. Reliability-based design optimization for crashworthiness of vehicle side impact. Structural and Multidisciplinary Optimization, 26:272–283, 2004. Samuel Daulton, Maximilian Balandat, and Eytan Bakshy. Parallel Bayesian optimization of multiple noisy objectives with expected hypervolume improvement. In Advances in Neural Information Processing Systems, volume 34. Curran Associates, Inc., 2021. Roman Garnett. Bayesian optimization. Cambridge University Press, 2023. Yaochu Jin. Surrogate-assisted evolutionary computation: Recent advances and future challenges. Swarm and Evolutionary Computation, 1(2):61–70, 2011. Jing Liang, Yahang Lou, Mingyuan Yu, Ying Bi, and Kunjie Yu. A survey of surrogate-assisted evolutionary algorithms for expensive optimization. Journal of Membrane Computing, pages 1–20, 2024. Joshua Knowles. Parego: A hybrid algorithm with on-line landscape approximation for expensive multiobjective optimization problems. IEEE transactions on Evolutionary Computation, 10(1):50–66, 2006. 10

arXiv Template

A P REPRINT

Biswajit Paria, Kirthevasan Kandasamy, and Barnabás Póczos. A flexible framework for multi-objective Bayesian optimization using random scalarizations. In Ryan P. Adams and Vibhav Gogate, editors, Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, volume 115 of Proceedings of Machine Learning Research, pages 766–776. PMLR, 2020. Xi Lin, Zhiyuan Yang, Xiaoyuan Zhang, and Qingfu Zhang. Pareto set learning for expensive multi-objective optimization. In Advances in Neural Information Processing Systems, volume 35, pages 19231–19247. Curran Associates, Inc., 2022. Samuel Daulton, Maximilian Balandat, and Eytan Bakshy. Differentiable expected hypervolume improvement for parallel multi-objective Bayesian optimization. In Advances in Neural Information Processing Systems, volume 33, pages 9851–9864. Curran Associates, Inc., 2020. Michael TM Emmerich, Kyriakos C Giannakoglou, and Boris Naujoks. Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels. IEEE Transactions on Evolutionary Computation, 10 (4):421–439, 2006. Ben Tu, Axel Gandy, Nikolas Kantas, and Behrang Shafei. Joint entropy search for multi-objective Bayesian optimization. In Advances in Neural Information Processing Systems, volume 35, pages 9922–9938. Curran Associates, Inc., 2022. Indraneel Das and John E Dennis. Normal-boundary intersection: A new method for generating the Pareto surface in nonlinear multicriteria optimization problems. SIAM journal on optimization, 8(3):631–657, 1998. Bingdong Li, Zixiang Di, Yongfan Lu, Hong Qian, Feng Wang, Peng Yang, Ke Tang, and Aimin Zhou. Expensive multi-objective Bayesian optimization based on diffusion models. Proceedings of the AAAI Conference on Artificial Intelligence, 39(25):27063–27071, 2025. Eckart Zitzler and Lothar Thiele. Multiobjective evolutionary algorithms: A comparative case study and the strength Pareto approach. IEEE Transactions on Evolutionary Computation, 3(4):257–271, 1999. Maximilian Balandat, Brian Karrer, Daniel Jiang, Samuel Daulton, Ben Letham, Andrew G Wilson, and Eytan Bakshy. Botorch: A framework for efficient monte-carlo Bayesian optimization. In Advances in Neural Information Processing Systems, volume 33, pages 21524–21538. Curran Associates, Inc., 2020. Mina Konakovic Lukovic, Yunsheng Tian, and Wojciech Matusik. Diversity-guided multi-objective Bayesian optimization with batch evaluations. In Advances in Neural Information Processing Systems, volume 33, pages 17708–17720. Curran Associates, Inc., 2020. Samuel Daulton, Sait Cakmak, Maximilian Balandat, Michael A Osborne, Enlu Zhou, and Eytan Bakshy. Robust multi-objective Bayesian optimization under input noise. In International Conference on Machine Learning, pages 4831–4866. PMLR, 2022a. Tinkle Chugh. Scalarizing functions in Bayesian multiobjective optimization. In 2020 IEEE Congress on Evolutionary Computation (CEC), pages 1–8. IEEE, 2020. Tze Leung Lai and Herbert Robbins. Asymptotically efficient adaptive allocation rules. Advances in Applied Mathematics, 6(1):4–22, 1985. Kaisa Miettinen. Nonlinear Multiobjective Optimization, volume 12. Springer Science & Business Media, 1999. Donald R Jones, Matthias Schonlau, and William J Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13(4):455–492, 1998. William R Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3/4):285–294, 1933. Ke Shang, Hisao Ishibuchi, Linjun He, and Lie Meng Pang. A survey on the hypervolume indicator in evolutionary multiobjective optimization. IEEE Transactions on Evolutionary Computation, 25(1):1–20, 2020. Ivo Couckuyt, Dirk Deschrijver, and Tom Dhaene. Fast calculation of multiobjective probability of improvement and expected improvement criteria for Pareto optimization. Journal of Global Optimization, 60(3):575–594, 2014. David Gaudrie, Rodolphe Le Riche, Victor Picheny, Benoit Enaux, and Vincent Herbert. Budgeted multi-objective optimization with a focus on the central part of the Pareto front–extended version. arXiv preprint arXiv:1809.10482, 2018. David Gaudrie, Rodolphe Le Riche, Victor Picheny, Benoit Enaux, and Vincent Herbert. Targeting solutions in Bayesian multi-objective optimization: sequential and batch versions. Annals of Mathematics and Artificial Intelligence, 88(1): 187–212, 2020. Mickael Binois, Victor Picheny, Patrick Taillandier, and Abderrahmane Habbal. The kalai-smorodinsky solution for many-objective Bayesian optimization. Journal of Machine Learning Research, 21(150):1–42, 2020. 11

arXiv Template

A P REPRINT

Majid Abdolshah, Alistair Shilton, Santu Rana, Sunil Gupta, and Svetha Venkatesh. Multi-objective Bayesian optimisation with preferences over objectives. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. Raul Astudillo and Peter Frazier. Multi-attribute Bayesian optimization with interactive preference learning. In International Conference on Artificial Intelligence and Statistics, pages 4496–4507. PMLR, 2020. Ryota Ozaki, Kazuki Ishikawa, Youhei Kanzaki, Shion Takeno, Ichiro Takeuchi, and Masayuki Karasuyama. Multiobjective Bayesian optimization with active preference learning. Proceedings of the AAAI conference on artificial intelligence, 38(13):14490–14498, 2024. Joshua Hang Sai Ip, Ankush Chakrabarty, Ali Mesbah, and Diego Romeres. User preference meets Pareto-optimality in multi-objective Bayesian optimization. Proceedings of the AAAI Conference on Artificial Intelligence, 39(19): 20246–20254, 2025. Jürgen Branke, Kalyanmoy Deb, Kaisa Miettinen, and Roman Słowiński, editors. Multiobjective Optimization: Interactive and Evolutionary Approaches, volume 5252. Springer, 2008. Jean-Pierre Imhof. Computing the distribution of quadratic forms in normal variables. Biometrika, 48(3/4):419–426, 1961. Harold Ruben. Probability content of regions under spherical normal distributions, iv: The distribution of homogeneous and non-homogeneous quadratic functions of normal variables. The Annals of Mathematical Statistics, 33(2): 542–570, 1962. Abhranil Das. New methods to compute the generalized chi-square distribution. Journal of Statistical Computation and Simulation, pages 1–35, 2025. Anton J Kleywegt, Alexander Shapiro, and Tito Homem-de Mello. The sample average approximation method for stochastic discrete optimization. SIAM Journal on Optimization, 12(2):479–502, 2002. Ilya M Sobol. The distribution of points in a cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics, 7(4):86–112, 1967. Sebastian Ament, Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy. Unexpected improvements to expected improvement for Bayesian optimization. In Advances in Neural Information Processing Systems, volume 36, pages 20577–20612. Curran Associates, Inc., 2023. Kalyanmoy Deb, Lothar Thiele, Marco Laumanns, and Eckart Zitzler. Scalable test problems for evolutionary multiobjective optimization. In Evolutionary multiobjective optimization: theoretical advances and applications, pages 105–145. Springer, 2005. Ran Cheng, Miqing Li, Ye Tian, Xingyi Zhang, Shengxiang Yang, Yaochu Jin, and Xin Yao. A benchmark test suite for evolutionary many-objective optimization. Complex & Intelligent Systems, 3:67–81, 2017. Miqing Li, Shengxiang Yang, and Xiaohui Liu. Shift-based density estimation for pareto-based algorithms in manyobjective optimization. IEEE Transactions on Evolutionary Computation, 18(3):348–365, 2014a. Ke Li, Kalyanmoy Deb, Qingfu Zhang, and Sam Kwong. An evolutionary many-objective optimization algorithm based on dominance and decomposition. IEEE transactions on Evolutionary Computation, 19(5):694–716, 2014b. Miqing Li, Shengxiang Yang, and Xiaohui Liu. Bi-goal evolution for many-objective optimization problems. Artificial Intelligence, 228:45–65, 2015b. Ryoji Tanabe and Hisao Ishibuchi. An easy-to-use real-world multi-objective optimization problem suite. Applied Soft Computing, 89:106078, 2020. Daniel Hernandez-Lobato, Jose Hernandez-Lobato, Amar Shah, and Ryan Adams. Predictive entropy search for multi-objective Bayesian optimization. In Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1492–1501, New York, USA, 2016. PMLR. Chao Jiang and Miqing Li. Trading off quality and uncertainty through multi-objective optimisation in batch Bayesian optimisation. Proceedings of the AAAI Conference on Artificial Intelligence, 39(25):27027–27035, 2025a. Miqing Li and Xin Yao. Quality evaluation of solution sets in multiobjective optimisation: A survey. ACM Computing Surveys (CSUR), 52(2):1–38, 2019. Miqing Li, Tao Chen, and Xin Yao. How to evaluate solutions in pareto-based search-based software engineering: A critical review and methodological guidance. IEEE Transactions on Software Engineering, 48(5):1771–1799, 2022. Hisao Ishibuchi, Ryo Imada, Yu Setoguchi, and Yusuke Nojima. How to specify a reference point in hypervolume calculation for fair performance comparison. Evolutionary Computation, 26(3):411–440, 2018. 12

arXiv Template

A P REPRINT

Frank Wilcoxon. Individual comparisons by ranking methods. In Breakthroughs in Statistics: Methodology and Distribution, pages 196–202. Springer, 1992. Sture Holm. A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, pages 65–70, 1979. Shinkyu Jeong and Shigeru Obayashi. Efficient global optimization (ego) for multi-objective problem and data mining. In 2005 IEEE congress on evolutionary computation, volume 3, pages 2138–2145. IEEE, 2005. Dianne Carrol Bautista. A sequential design for approximating the Pareto front using the expected Pareto improvement function. The Ohio State University, 2009. Joshua Svenson. Computer experiments: Multiobjective optimization and sensitivity analysis. PhD thesis, The Ohio State University, 2011. Joshua Svenson and Thomas Santner. Multiobjective optimization of expensive-to-evaluate deterministic computer simulator models. Computational Statistics & Data Analysis, 94:250–264, 2016. Marcela Zuluaga, Andreas Krause, and Markus Püschel. ϵ-pal: An active learning approach to the multi-objective optimization problem. Journal of Machine Learning Research, 17(104):1–32, 2016. URL http://jmlr.org/ papers/v17/15-047.html. Dawei Zhan, Yuansheng Cheng, and Jun Liu. Expected improvement matrix-based infill criteria for expensive multiobjective optimization. IEEE Transactions on Evolutionary Computation, 21(6):956–975, 2017. Victor Picheny, Mickael Binois, and Abderrahmane Habbal. A Bayesian optimization approach to find nash equilibria. Journal of Global Optimization, 73:171–192, 2019. Syrine Belakaria, Aryan Deshwal, Nitthilan Kannappan Jayakodi, and Janardhan Rao Doppa. Uncertainty-aware search framework for multi-objective Bayesian optimization. Proceedings of the AAAI Conference on Artificial Intelligence, 34(06):10044–10052, 2020a. Gustavo Malkomes, Bolong Cheng, Eric H Lee, and Mike Mccourt. Beyond the Pareto efficient frontier: Constraint active search for multiobjective experimental design. In Proceedings of the 38th International Conference on Machine Learning, volume 139, pages 7423–7434. PMLR, 2021. Ji Won Park, Nataša Tagasovska, Michael Maser, Stephen Ra, and Kyunghyun Cho. Botied: Multi-objective bayesian optimization with tied multivariate ranks. arXiv preprint arXiv:2306.00344, 2023. Lam Ngo, Huong Ha, Jeffrey Chan, and Hongyu Zhang. Mobo-osd: Batch multi-objective Bayesian optimization via orthogonal search directions. arXiv preprint arXiv:2510.20872, 2025. Andy J Keane. Statistical improvement criteria for use in multiobjective design optimization. AIAA journal, 44(4): 879–891, 2006. James Parr. Improvement criteria for constraint handling and multiobjective optimization. PhD thesis, University of Southampton, 2013. Marcela Zuluaga, Guillaume Sergent, Andreas Krause, and Markus Püschel. Active learning for multi-objective optimization. In International conference on machine learning, pages 462–470. PMLR, 2013. Alaleh Ahmadianshalchi, Syrine Belakaria, and Janardhan Rao Doppa. Pareto front-diverse batch multi-objective Bayesian optimization. Proceedings of the AAAI Conference on Artificial Intelligence, 38(10):10784–10794, 2024. Liang Zhao and Qingfu Zhang. Exact formulas for the computation of expected tchebycheff improvement. In 2023 IEEE Congress on Evolutionary Computation (CEC), pages 1–8. IEEE, 2023. Qingfu Zhang, Wudong Liu, Edward Tsang, and Botond Virginas. Expensive multiobjective optimization by MOEA/D with Gaussian process model. IEEE Transactions on Evolutionary Computation, 14(3):456–474, 2010. Nobuo Namura, Koji Shimoyama, and Shigeru Obayashi. Expected improvement of penalty-based boundary intersection for expensive multiobjective optimization. IEEE Transactions on Evolutionary Computation, 21(6):898–913, 2017. Richard Zhang and Daniel Golovin. Random hypervolume scalarizations for provable multi-objective black box optimization. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 11096–11105. PMLR, 2020. Diantong Li, Fengxue Zhang, Chong Liu, and Yuxin Chen. Constrained multi-objective Bayesian optimization through optimistic constraints estimation. arXiv preprint arXiv:2411.03641, 2024. Wolfgang Ponweiser, Tobias Wagner, Dirk Biermann, and Markus Vincze. Multiobjective optimization on a limited budget of evaluations using model-assisted S-metric selection. In Parallel Problem Solving from Nature – PPSN X, pages 784–794. Springer, 2008. 13

arXiv Template

A P REPRINT

Ashwin Renganathan and Kade Carlson. qPOTS: Efficient batch multiobjective Bayesian optimization via pareto optimal Thompson sampling. In Proceedings of The 28th International Conference on Artificial Intelligence and Statistics, volume 258 of Proceedings of Machine Learning Research, pages 4051–4059. PMLR, 03–05 May 2025. Shriya Bhatija, Paul-David Zuercher, Jakob Thumm, and Thomas Bohné. Multi-objective causal bayesian optimization. arXiv preprint arXiv:2502.14755, 2025. Sam Daulton, Maximilian Balandat, and Eytan Bakshy. Hypervolume knowledge gradient: A lookahead approach for multi-objective Bayesian optimization with partial information. In Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 7167–7204. PMLR, 2023. Jixiang Qing, Ivo Couckuyt, and Tom Dhaene. A robust multi-objective Bayesian optimization framework considering input uncertainty. Journal of Global Optimization, 86(3):693–711, 2023. Kaifeng Yang, Michael Emmerich, André Deutz, and Thomas Bäck. Multi-objective Bayesian global optimization using expected hypervolume improvement gradient. Swarm and Evolutionary Computation, 44:945–956, 2019a. Kaifeng Yang, Michael Emmerich, André Deutz, and Thomas Bäck. Efficient computation of expected hypervolume improvement using box decomposition algorithms. Journal of Global Optimization, 75(1):3–34, 2019b. Jingda Deng, Jianyong Sun, Qingfu Zhang, and Hui Li. Expected hypervolume improvement is a particular hypervolume improvement. Proceedings of the AAAI Conference on Artificial Intelligence, 39(15):16217–16225, 2025. Eduardo C Garrido-Merchán, Daniel Fernández-Sánchez, and Daniel Hernández-Lobato. Parallel predictive entropy search for multi-objective Bayesian optimization with constraints applied to the tuning of machine learning algorithms. Expert Systems with Applications, 215:119328, 2023. Eduardo C. Garrido-Merchán and Daniel Hernández-Lobato. Predictive entropy search for multi-objective Bayesian optimization with constraints. Neurocomputing, 361:50–68, 2019. Syrine Belakaria, Aryan Deshwal, and Janardhan Rao Doppa. Max-value entropy search for multi-objective Bayesian optimization. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. Syrine Belakaria, Aryan Deshwal, and Janardhan Rao Doppa. Output space entropy search framework for multiobjective Bayesian optimization. Journal of artificial intelligence research, 72:667–715, 2021. Shinya Suzuki, Shion Takeno, Tomoyuki Tamura, Kazuki Shitara, and Masayuki Karasuyama. Multi-objective Bayesian optimization using Pareto-frontier entropy. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 9279–9288. PMLR, 2020. Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrixmatrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018. Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019. Himanshu Jain and Kalyanmoy Deb. An improved adaptive approach for elitist nondominated sorting genetic algorithm for many-objective optimization. In International Conference on Evolutionary Multi-Criterion Optimization, pages 307–321. Springer, 2013b. Chao Jiang and Miqing Li. Multi-objectivising acquisition functions in Bayesian optimisation. ACM Transactions on Evolutionary Learning and Optimization, 5(2), 2025b. Miguel González-Duque, Richard Michael, Simon Bartels, Yevgen Zainchkovskyy, Sø ren Hauberg, and Wouter Boomsma. A survey and benchmark of high-dimensional Bayesian optimization of discrete sequences. In Advances in Neural Information Processing Systems, volume 37, pages 140478–140508. Curran Associates, Inc., 2024. Leonard Papenmeier, Luigi Nardi, and Matthias Poloczek. Bounce: Reliable high-dimensional Bayesian optimization for combinatorial and mixed spaces. In Advances in Neural Information Processing Systems, volume 36, pages 1764–1793. Curran Associates, Inc., 2023. Maria Laura Santoni, Elena Raponi, Renato De Leone, and Carola Doerr. Comparison of high-dimensional Bayesian optimization algorithms on BBOB. ACM Transactions on Evolutionary Learning and Optimisation, 4(3):1–33, 2024. Ziyu Wang, Frank Hutter, Masrour Zoghi, David Matheson, and Nando De Feitas. Bayesian optimization in a billion dimensions via random embeddings. Journal of Artificial Intelligence Research, 55:361–387, 2016. Vu Viet Hoang, Hung The Tran, Sunil Gupta, and Vu Nguyen. High dimensional Bayesian optimization using lasso variable selection. arXiv preprint arXiv:2504.01743, 2025. 14

arXiv Template

A P REPRINT

Mickael Binois and Nathan Wycoff. A survey on high-dimensional Gaussian process modeling with application to Bayesian optimization. ACM Transactions on Evolutionary Learning and Optimization, 2(2):1–26, 2022. Taicai Chen, Yue Duan, Dong Li, Lei Qi, Yinghuan Shi, and Yang Gao. PG-LBO: enhancing high-dimensional Bayesian optimization with pseudo-label and Gaussian process guidance. Proceedings of the AAAI Conference on Artificial Intelligence, 38(10):11381–11389, 2024. Amin Nayebi, Alexander Munteanu, and Matthias Poloczek. A framework for Bayesian optimization in embedded subspaces. In International Conference on Machine Learning, pages 4752–4761. PMLR, 2019. Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka. Batched large-scale Bayesian optimization in high-dimensional spaces. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84, pages 745–754. PMLR, 2018. Zhitong Xu, Haitao Wang, Jeff M Phillips, and Shandian Zhe. Standard Gaussian process is all you need for highdimensional Bayesian optimization. In The Thirteenth International Conference on Learning Representations, 2025. David Eriksson and Martin Jankowiak. High-dimensional Bayesian optimization with sparse axis-aligned subspaces. In Uncertainty in Artificial Intelligence, pages 493–503. PMLR, 2021. Ian Delbridge, David Bindel, and Andrew Gordon Wilson. Randomly projected additive Gaussian processes for regression. In International Conference on Machine Learning, pages 2453–2463. PMLR, 2020. Eric Han, Ishank Arora, and Jonathan Scarlett. High-dimensional Bayesian optimization via tree-structured additive models. Proceedings of the AAAI Conference on Artificial Intelligence, 35(9):7630–7638, 2021. Juliusz Krzysztof Ziomek and Haitham Bou Ammar. Are random decompositions all we need in high dimensional Bayesian optimisation? In International Conference on Machine Learning, pages 43347–43368. PMLR, 2023. Kirill Antonov, Elena Raponi, Hao Wang, and Carola Doerr. High dimensional Bayesian optimization with kernel principal component analysis. In International Conference on Parallel Problem Solving from Nature, pages 118–131. Springer, 2022. Ben Letham, Roberto Calandra, Akshara Rai, and Eytan Bakshy. Re-examining linear embeddings for high-dimensional Bayesian optimization. In Advances in Neural Information Processing Systems, volume 33, pages 1546–1558. Curran Associates, Inc., 2020. Elena Raponi, Hao Wang, Mariusz Bujny, Simonetta Boria, and Carola Doerr. High dimensional Bayesian optimization assisted by principal component analysis. In Parallel Problem Solving from Nature–PPSN XVI: 16th International Conference, PPSN 2020, Leiden, The Netherlands, September 5-9, 2020, Proceedings, Part I 16, pages 169–183. Springer, 2020. Samuel Daulton, David Eriksson, Maximilian Balandat, and Eytan Bakshy. Multi-objective Bayesian optimization over high-dimensional search spaces. In Uncertainty in Artificial Intelligence, pages 507–517. PMLR, 2022b. Youssef Diouane, Victor Picheny, Rodolophe Le Riche, and Alexandre Scotto Di Perrotolo. Trego: a trust-region framework for efficient global optimization. Journal of Global Optimization, 86(1):1–23, 2023. David Eriksson, Michael Pearce, Jacob Gardner, Ryan D Turner, and Matthias Poloczek. Scalable global optimization via local Bayesian optimization. In Advances in Neural Information Processing Systems, volume 32, pages 5496–5507. Curran Associates, Inc., 2019. Carl Hvarfner, Erik Orm Hellsten, and Luigi Nardi. Vanilla Bayesian optimization performs great in high dimensions. In Proceedings of the 41st International Conference on Machine Learning, volume 235, pages 20793–20817. PMLR, 21–27 Jul 2024. Leonard Papenmeier, Matthias Poloczek, and Luigi Nardi. Understanding high-dimensional Bayesian optimization. arXiv preprint arXiv:2502.09198, 2025. Syrine Belakaria, Aryan Deshwal, and Janardhan Rao Doppa. Multi-fidelity multi-objective Bayesian optimization: An output space entropy search approach. Proceedings of the AAAI Conference on artificial intelligence, 34(06): 10035–10043, 2020b. Kirthevasan Kandasamy, Gautam Dasarathy, Jeff Schneider, and Barnabás Póczos. Multi-fidelity Bayesian optimisation with continuous approximations. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 1799–1808. PMLR, 2017. Shibo Li, Wei Xing, Robert Kirby, and Shandian Zhe. Multi-fidelity Bayesian optimization via deep neural networks. In Advances in Neural Information Processing Systems, volume 33, pages 8521–8531. Curran Associates, Inc., 2020. Henry B. Moss, David S. Leslie, Javier Gonzalez, and Paul Rayson. Gibbon: General-purpose information-based Bayesian optimisation. Journal of Machine Learning Research, 22(235):1–49, 2021. 15

arXiv Template

A P REPRINT

Jialin Song, Yuxin Chen, and Yisong Yue. A general framework for multi-fidelity Bayesian optimization with Gaussian processes. In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, volume 89 of Proceedings of Machine Learning Research, pages 3158–3167. PMLR, 2019. Shion Takeno, Hitoshi Fukuoka, Yuhki Tsukada, Toshiyuki Koyama, Motoki Shiga, Ichiro Takeuchi, and Masayuki Karasuyama. Multi-fidelity Bayesian optimization with max-value entropy search and its parallelization. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 9334–9345. PMLR, 2020. Jian Wu, Saul Toscano-Palmerin, Peter I. Frazier, and Andrew Gordon Wilson. Practical multi-fidelity Bayesian optimization for hyperparameter tuning. In Proceedings of The 35th Uncertainty in Artificial Intelligence Conference, volume 115 of Proceedings of Machine Learning Research, pages 788–798. PMLR, 22–25 Jul 2020. Yehong Zhang, Trong Nghia Hoang, Bryan Kian Hsiang Low, and Mohan Kankanhalli. Information-based multi-fidelity Bayesian optimization. In NIPS workshop on Bayesian optimization, volume 49. Journal of Machine Learning Research JMLR. org Cambridge, MA, 2017.

arXiv Template

A P REPRINT

Appendix to: Do We Really Need to Approach the Entire Pareto Front in Many-Objective Bayesian Optimisation? A

Multi-objective Bayesian Optimisation (MOBO)

MOBO consists of two main steps, i.e., training m Gaussian process models based on the observed solutions and optimising an acquisition function α(x) : X → R to select a solution for evaluation. In this work, we model each objective with an independent Gaussian process fi ∼ GP(mi (x), ki (x, x′ )), where mi (x) : X → R is the ith mean function, and ki (·, ·) : X × X → R is the ith covariance function. We use the notation K(A, B) to represent the covariance matrix at all pairs of solutions in set A and in set B. Given n observed solutions Dn = {(xt , y t )}nt=1 where y t = f (xt ) + ζ t and the noise ζ t ∼ N (0, diag(σζ2 )), the posterior distribution of ith objective at a new location x is a Gaussian distribution: p(fi (x)|Dn ) ∼ N (µi (x), σi2 (x)) (6) µi (x) = K(x, X n )(K(X n , X n ) + σζ2i I)−1 Yi

(7)

σi2 (x) = K(x, x) − K(x, X n )((K(X n , X n ) + σζ2i I)−1 K(X n , x)

(8)

where µi (x) and σi2 (x) are the mean and variance at x, respectively; X n = (x1 , . . . , xn ) ∈ Rn×d and Yin = (yi1 , . . . , yin ) ∈ Rn are the matrix of evaluated solutions and the corresponding vector of y values, respectively; σζ2i is the variance of the observation noise ζi ∼ N (0, σζ2i ), and corresponds to the i-th diagonal entry of the noise covariance matrix diag(σζ2 );

Illustrative Example of ESPI

To help understand the proposed ESPI, an illustration of ESPI in a bi-objective case is shown in Figure 5. In this example, the utopian point is denoted by the black star, while the red dots represent the nondominated solutions from the current dataset. A new candidate point, whose objective values are yet to be observed, is shown as the blue dot. The red dotted line indicates the shortest Euclidean distance g ∗ from the utopian point to the existing nondominated set. In contrast, the blue dotted line represents the distance from the new point to the utopian point. When a new point lies within the shaded blue region and is closer to the utopian point, the improvement ISP (·) is higher. Nondominated points in n New point (not observed) Utopian point

p(f2(x)|Dn)

Figure 5: Illustration of the proposed expected single-point improvement (ESPI) in a bi-objective space. The figure shows the utopian point (black star), nondominated points in the current dataset Dn (red dots), and a new candidate point whose true objective values are not yet observed (blue dot). The dashed arc represents the current best Euclidean distance from the utopian point to the nondominated set, denoted as g ∗ (red dotted line). The blue dotted line represents the distance from the new candidate point to the utopian point z ∗ , which may improve upon the current best distance g ∗ . When a new point lies within the shaded blue region and is closer to the utopian point, the improvement ISP (·) is higher.

ISP (·) > 0

g∗ p(f1(x)|Dn)

arXiv Template

A P REPRINT

Extended Related Work

Over the last two decades, a variety of BO methods have been proposed to tackle expensive multi-objective optimisation problems [Jeong and Obayashi, 2005, Bautista, 2009, Svenson, 2011, Svenson and Santner, 2016, Zuluaga et al., 2016, Zhan et al., 2017, Picheny et al., 2019, Belakaria et al., 2020a, Malkomes et al., 2021, Park et al., 2023, Ngo et al., 2025]. Most of them aim to identify a good approximation of the entire Pareto front [Keane, 2006, Svenson, 2011, Parr, 2013, Zuluaga et al., 2013, Ahmadianshalchi et al., 2024]. To do so, some studies convert a multi-objective problem into multiple single-objective problems by a scalarisation function (e.g., random augmented Tchebycheff scalarisation). They then optimise acquisition functions from single-objective BO to determine the next evaluation point(s) [Knowles, 2006, Zhao and Zhang, 2023]. In these studies, different acquisition functions are employed, such as expected improvement (EI) [Jones et al., 1998] in Knowles [2006], Zhang et al. [2010], Namura et al. [2017], Chugh [2020], Thompson sampling (TS) [Thompson, 1933] in Paria et al. [2020], Zhang and Golovin [2020], and upper confidence bound (UCB) [Lai and Robbins, 1985] in Paria et al. [2020], Zhang and Golovin [2020], Li et al. [2024]. The remaining studies directly optimise multi-objective problems by considering the definition of optimality in multiobjective optimisation, i.e., the Pareto dominance relation. A representative approach is to use HV since maximising the HV value is equivalent to finding the entire Pareto front [Ponweiser et al., 2008, Couckuyt et al., 2014, Daulton et al., 2020, 2021, Renganathan and Carlson, 2025, Bhatija et al., 2025]. Along this line, expected hypervolume improvement (EHVI) is widely considered [Emmerich et al., 2006, Daulton et al., 2023, Qing et al., 2023, Yang et al., 2019a,b, Deng et al., 2025] as it is a natural extension of the EI for multi-objective optimisation. Another idea is to leverage information theory to guide exploration toward regions likely contributing to the Pareto front. Such methods focus on improving the posterior of optimal inputs (i.e., the approximated Pareto set) [Garrido-Merchán et al., 2023, Garrido-Merchán and Hernández-Lobato, 2019, Hernandez-Lobato et al., 2016], optimal outputs (i.e., the approximated Pareto front) [Belakaria et al., 2019, 2021, Suzuki et al., 2020], or both of them [Tu et al., 2022]. That said, there do exist a few studies that do not aim to identify the entire Pareto front. Among them, some attempt to use decision-maker preferences to guide the search towards specific region(s) [Abdolshah et al., 2019, Astudillo and Frazier, 2020, Ozaki et al., 2024, Ip et al., 2025]. Such methods typically adjust the target region by eliciting or updating decision-maker preferences during optimisation. Another attempt is to directly target certain region(s) of the Pareto front, without an assumption that decision-maker preferences can be available or elicited [Gaudrie et al., 2018, 2020, Binois et al., 2020]. For example, Gaudrie et al. [2018, 2020] propose the Centred Expected Hypervolume Improvement (C-EHVI) by dynamically adjusting the reference point using the Kalai-Smorodinsky equilibrium (also used in Binois et al. [2020]) to approach the central part of the Pareto front. Such work is more relevant to our study, and we have thus included C-EHVI in our experimental comparison.

Theoretical Results

D.1

Proof of Theorem 4.1.

We consider the setting from Balandat et al. [2020, Section D.5]. Let ϵt ∼ N (0, Im ).4 Using the reparameterisation trick, we can write the posterior at x as f¯t (x, ϵt ) = µ̄(x) + L̄(x)ϵt where µ̄(x) : Rd → Rm is the multi-output GP’s posterior mean; L̄(x) ∈ Rm×m is a root decomposition (often a Cholesky decomposition) of the multi-output GP’s posterior covariance K̄ ∈ Rm×m ; and ϵt ∈ Rm . Let Ā(x, ϵt ) = max 0, g ∗ − ∥f¯t (x, ϵt ) − z ∗ ∥ ∗ where g ∗ = minx∈X̄ n g(f (x), z ∗ ) and z ∗ = (z1∗ , z2∗ , . . . , zm ) is the utopian point. Following Balandat et al. [2020, Theorem 3], we need to show that there exists an integrable function ℓ : Rm 7→ R such that for almost every ϵt and all x, y ∈ X ⊂ Rd ,

|Ā(x, ϵt ) − Ā(y, ϵt )| ≤ ℓ(ϵt )∥x − y∥. 4

(9)

Theorem 4.1 can be extended to handle non-iid base samples from a family of quasi-Monte Carlo methods as in Balandat et al. [2020].

arXiv Template

A P REPRINT

Let Ā(x, ϵt ) = max 0, g ∗ − ∥f¯t (x, ϵt ) − z ∗ ∥ 1 ∗ g − ∥f¯t (x, ϵt ) − z ∗ ∥ + g ∗ − ∥f¯t (x, ϵt ) − z ∗ ∥ . = 2 Hence we obtain, Ā(x, ϵt ) − Ā(y, ϵt ) 1 1 ¯ ∥ft (y, ϵt ) − z ∗ ∥ − ∥f¯t (x, ϵt ) − z ∗ ∥ + g ∗ − ∥f¯t (x, ϵt ) − z ∗ ∥ − g ∗ − ∥f¯t (y, ϵt ) − z ∗ ∥ . = 2 2 Let I1 := ∥f¯t (y, ϵt ) − z ∗ ∥ − ∥f¯t (x, ϵt ) − z ∗ ∥ and I2 := g ∗ − ∥f¯t (x, ϵt ) − z ∗ ∥ − g ∗ − ∥f¯t (y, ϵt ) − z ∗ ∥ . Hence Ā(x, ϵt ) − Ā(y, ϵt ) ≤

1 1 |I1 | + |I2 |. 2 2

We observe that |I1 | = ∥f¯t (y, ϵt ) − z ∗ ∥ − ∥f¯t (x, ϵt ) − z ∗ ∥ ≤ ∥f¯t (y, ϵt ) − f¯t (x, ϵt )∥ = ∥µ̄(y) + L̄(y)ϵt − (µ̄(x) + L̄(x)ϵt )∥

≤ ∥µ̄(y) − µ̄(x)∥ + ∥(L̄(y) − L̄(x))ϵt ∥. Since X is compact, and µ̄ and L̄ have uniformly bounded gradients, they are Lipschitz. There exist Cµ1 , CL1 < ∞ such that |I1 | ≤ ∥µ̄(y) − µ̄(x)∥ + ∥(L̄(y) − L̄(x))ϵt ∥ ≤ ℓI1 (ϵt )∥x − y∥

where ℓI1 (ϵt ) := Cµ1 + CL1 ∥ϵt ∥. Furthermore, |I2 | = g ∗ − ∥f¯t (x, ϵt ) − z ∗ ∥ − g ∗ − ∥f¯t (y, ϵt ) − z ∗ ∥ ≤ ∥f¯t (x, ϵt ) − z ∗ ∥ − ∥f¯t (y, ϵt ) − z ∗ ∥ ≤ ∥f¯t (x, ϵt ) − f¯t (y, ϵt )∥ = ∥µ̄(x) + L̄(x)ϵt − (µ̄(y) + L̄(y)ϵt )∥

≤ ∥µ̄(x) − µ̄(y)∥ + ∥(L̄(x) − L̄(y))ϵt ∥.

Since X is compact, and µ̄ and L̄ have uniformly bounded gradients, they are Lipschitz. There exist Cµ2 , CL2 < ∞ such that |I2 | ≤ ℓI2 (ϵt )∥x − y∥ where ℓI2 (ϵt ) := Cµ2 + CL2 ∥ϵt ∥. Hence 1 1 |I1 | + |I2 | 2 2 1 1 ≤ ℓI1 (ϵt )∥x − y∥ + ℓI2 (ϵt )∥x − y∥ 2 2 1 = (ℓI1 (ϵt ) + ℓI2 (ϵt ))∥x − y∥. 2

Ā(x, ϵt ) − Ā(y, ϵt ) ≤

Hence Ā(x, ϵt ) − Ā(y, ϵt ) ≤ ℓ(ϵt )∥x − y∥ where ℓ(ϵt ) := (Cµ1 + Cµ2 ) + (CL1 + CL2 )∥ϵt ∥. Note that ℓ(ϵt ) is integrable because all absolute moments exist for the Gaussian distribution. Since this satisfies the criteria for Theorem 3 in Balandat et al. [2020], the theorem holds for ESPI. 19

arXiv Template

D.2

A P REPRINT

Theorem D.1 and Its Proof

Theorem D.1. Suppose that X is compact and that f has a multi-output GP prior with continuously differentiable mean and covariance functions. Let X n = {xt }nt=1 be the set of observed decision vectors in Dn , α∗NESPI := maxx∈X αNESPI (x) denote the maximum of NESPI, S ∗ := arg maxx∈X αNESPI (x) denote the set of maximisers t N of αNESPI , α̂N NESPI (x) denote the deterministic acquisition function via the base samples {ϵ }t=1 ∼ N (0, I(n+1)m ). ∗ N Suppose that x̂N ∈ arg maxx∈X α̂NESPI (x), then ∗ ∗ (1) α̂N NESPI (x̂N ) → αNESPI a.s.,

(2) d(x̂N∗ , S ∗ ) → 0 a.s., where d(x̂N∗ , S ∗ ) := inf x∈S ∗ ∥x̂N∗ − x∥. Proof of Theorem D.1. Let X n := [(x1 )T , . . . , (xn )T ]T ∈ Rnd be the n observed points; xn+1 ∈ X ⊂ Rd be a candidate point; and ϵt ∈ R(n+1)m with ϵt ∈ N (0, I(n+1)m ). Let f˜t (X n , xn+1 ) := [f˜t (x1 ), . . . , f˜t (xn ), f˜t (xn+1 )] denote the tth sample of the corresponding objectives, so that we can write the posterior via the reparameterisation trick as: ft (X n , xn+1 , ϵt ) = µ(X n , xn+1 ) + L(X n , xn+1 )ϵt where µ(X n , xn+1 ) : R(n+1)d → R(n+1)m is the multi-output GP’s posterior mean; L(X n , xn+1 ) ∈ R(n+1)m×(n+1)m is a root decomposition of the multi-output GP’s posterior covariance K ∈ R(n+1)m×(n+1)m . Let f (m) (xi , ϵt ) := Si µ(X n , xn+1 ) + L(X n , xn+1 )ϵt , where f (m) (xi , ϵt ) ∈ Rm represents the posterior at point xi and Si ∈ Rm×(n+1)m , i = 1, . . . , n + 1 is the selector matrix used to extract the corresponding element for the xi . Let A(xn+1 , ϵt ; X n ) = max 0, gˆt∗ (X n ) − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ . (m)

Let gˆt∗ (X n ) := minxobs ∈{xi }ni=1 ∥ft

(xobs , ϵt ) − z ∗ ∥, hence we obtain:

A(xn+1 , ϵt ; X n ) = max 0,

min n ∥f (m) (xobs , ϵt ) − z ∗ ∥ − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ .

xobs ∈{xi }i=1

Following Balandat et al. [2020, Theorem 3], we need to show that there exists an integrable function ℓ : Rm 7→ R such that for almost every ϵt and all xn+1 , y n+1 ∈ X ⊂ Rd , |A(xn+1 , ϵt ; X n ) − A(y n+1 , ϵt ; X n )| ≤ ℓ(ϵt )∥xn+1 − y n+1 ∥.

(10)

Let A(xn+1 , ϵt ; X n ) = max 0, gˆt∗ (X n ) − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ 1 ˆ∗ n = gt (X ) − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ + gˆt∗ (X n ) − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ . 2 Hence we obtain A(xn+1 , ϵt ; X n ) − A(y n+1 , ϵt ; X n ) 1 (m) n+1 t = ∥f (y , ϵ ) − z ∗ ∥ − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ 2 1 ˆ∗ n + gt (X ) − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ − gˆt∗ (X n ) − ∥f (m) (y n+1 , ϵt ) − z ∗ ∥ . 2 Let I1′ := ∥f (m) (y n+1 , ϵt )−z ∗ ∥−∥f (m) (xn+1 , ϵt )−z ∗ ∥ and I2′ := gˆt∗ (X n )−∥f (m) (xn+1 , ϵt )−z ∗ ∥ − gˆt∗ (X n )− ∥f (m) (y n+1 , ϵt ) − z ∗ ∥ . Hence A(xn+1 , ϵt ; X n ) − A(y n+1 , ϵt ; X n ) ≤ 20

1 ′ 1 |I1 | + |I2′ |. 2 2

arXiv Template

A P REPRINT

We observe that |I1′ | = ∥f (m) (y n+1 , ϵt ) − z ∗ ∥ − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ ≤ ∥f (m) (y n+1 , ϵt ) − f (m) (xn+1 , ϵt )∥

= ∥µ(m) (y n+1 ) + L(m) (y n+1 )ϵt − (µ(m) (xn+1 ) + L(m) (xn+1 )ϵt )∥

≤ ∥µ(m) (y n+1 ) − µ(m) (xn+1 )∥ + ∥(L(m) (y n+1 ) − L(m) (xn+1 ))ϵt ∥ Since X is compact, and µ(m) and L(m) have uniformly bounded gradients, they are Lipschitz. There exist Cµ′ 1 , CL′ 1 < ∞ such that |I1′ | ≤ ∥µ(m) (y n+1 ) − µ(m) (xn+1 )∥ + ∥(L(m) (y n+1 ) − L(m) (xn+1 ))ϵt ∥ ≤ ℓI1′ (ϵt )∥xn+1 − y n+1 ∥.

where ℓI1′ (ϵt ) := Cµ′ 1 + CL′ 1 ∥ϵt ∥. Furthermore,

|I2′ | = gˆt∗ (X n ) − ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ − gˆt∗ (X n ) − ∥f (m) (y n+1 , ϵt ) − z ∗ ∥ ≤ ∥f (m) (xn+1 , ϵt ) − z ∗ ∥ − ∥f (m) (y n+1 , ϵt ) − z ∗ ∥ ≤ ∥f (m) (xn+1 , ϵt ) − f (m) (y n+1 , ϵt )∥

= ∥µ(m) (xn+1 ) + L(m) (xn+1 )ϵt − (µ(m) (y n+1 ) + L(m) (y n+1 )ϵt )∥

≤ ∥µ(m) (xn+1 ) − µ(m) (y n+1 )∥ + ∥(L(m) (xn+1 ) − L(m) (y n+1 ))ϵt ∥

Since X is compact, and µ(m) and L(m) have uniformly bounded gradients, they are Lipschitz. There exist Cµ′ 2 , CL′ 2 < ∞ such that where ℓI2′ (ϵt ) := Cµ′ 2 + CL′ 2 ∥ϵt ∥. Hence

|I2′ | ≤ ℓI2′ (ϵt )∥xn+1 − y n+1 ∥

1 ′ 1 |I | + |I2′ | 2 1 2 1 1 ≤ ℓI1′ (ϵt )∥xn+1 − y n+1 ∥ + ℓI2′ (ϵt )∥xn+1 − y n+1 ∥ 2 2 1 = (ℓI1′ (ϵt ) + ℓI2′ (ϵt ))∥xn+1 − y n+1 ∥ 2

A(xn+1 , ϵt ) − A(y n+1 , ϵt ) ≤

Hence A(xn+1 , ϵt ) − A(y n+1 , ϵt ) ≤ ℓ(ϵt )∥xn+1 − y n+1 ∥

where ℓ(ϵt ) := (Cµ′ 1 + Cµ′ 2 ) + (CL′ 1 + CL′ 2 )∥ϵt ∥. Note that ℓ(ϵt ) is integrable because all absolute moments exist for the Gaussian distribution. Since this satisfies the criteria for Theorem 3 in Balandat et al. [2020], the theorem holds for NESPI. It is worth noting that Theorems 4.1 and D.1 readily extend to the batch setting. It can be shown that the gradient of α̂ESPI (x) or α̂NESPI (x) is an unbiased estimator of the true gradient of αESPI or αNESPI , though it is not necessary for the SAA approach [Daulton et al., 2021].

Experiment Settings

E.1

Implementation details

All experiments were conducted using Python 3.12, with all methods developed on the open-source Python framework BoTorch [Balandat et al., 2020], which builds on GPyTorch [Gardner et al., 2018] for Gaussian process modelling and PyTorch [Paszke et al., 2019] for automatic differentiation. The computational studies were performed on a Red Hat Enterprise Linux 8.8 system, operating on a 64-bit x86 CPU architecture. The computing cluster utilised Intel Xeon Platinum 8360Y processors running at 2.40 GHz. The code for ParEGO, NParEGO, TS-TCH, EHVI, NEHVI, and JES is available at https://github.com/pytorch/ botorch. Specifically, we implement ParEGO and NParEGO based on the implementation provided by BoTorch.5 We 5

EI, NEI, and the BoTorch multi-objective tutorial.

arXiv Template

A P REPRINT

implement TS-TCH based on the implementation provided by BoTorch.6 We implement EHVI and NEHVI based on the implementation provided by BoTorch.7 We implement JES based on the implementation provided by BoTorch.8 The data and code are available at an anonymised repository for reproducibility: https://anonymous.4open. science/r/SPMO-B9DE. E.2

Method Details

For all the methods, we set the same 2(d + 1) points from a scrambled Sobol sequence and allow a maximum of 200 evaluations, following the practice in Daulton et al. [2020, 2021], Konakovic Lukovic et al. [2020]. All the methods use N = 128 Monte Carlo samples. For ParEGO [Knowles, 2006] and its noisy variant NParEGO [Daulton et al., 2021], we employ random scalarisations, whereby a weight vector w ∈ Rm is generated from P the unit simplex. The augmented Tchebycheff scalarisation function is applied, defined as g(y) = maxi (wi yi ) + α i (wi yi ). Log expected improvement and noisy log expected improvement are used as the acquisition functions in ParEGO and NParEGO, respectively, as recommended in Ament et al. [2023]. In the batch setting, q distinct weight vectors are sampled, and the acquisition function is optimised sequentially for each. Table 4: Reference points of the six benchmark problems used for hypervolume computation in EHVI and NEHVI, as well as for performance evaluation. Note that the referents points of the two real-world problems, i.e., car side impact design and car cab design, are set to (1.1, ..., 1.1) ∈ Rm in the normalised objective space following the practice in Tanabe and Ishibuchi [2020] the utopian and nadir points are available at https://github.com/ryojitanabe/ reproblems/tree/master/ideal_nadir_points). Problem

Reference point

Suggested by

DTLZ1 DTLZ2 Inverted DTLZ1 Inverted DTLZ2 Convex DTLZ2 Scaled DTLZ2 DTLZ3 DTLZ4 DTLZ5 DTLZ6 DTLZ7

(400.0, ..., 400.0) ∈ Rm (1.1, . . . , 1.1) ∈ Rm (400.0, . . . , 400.0) ∈ Rm (1.1, . . . , 1.1) ∈ Rm (1.1, . . . , 1.1) ∈ Rm (1.1 ∗ 20 , . . . , 1.1 ∗ 2m ) ∈ Rm (10000.0, . . . , 10000.0) ∈ Rm (1.1, . . . , 1.1) ∈ Rm (10.0, . . . , 10.0) ∈ Rm (10.0, . . . , 10.0) ∈ Rm (15.0, . . . , 15.0) ∈ Rm

Balandat et al. [2020], Chugh [2020] Balandat et al. [2020], Daulton et al. [2020], Ishibuchi et al. [2018] Chugh [2020], Ishibuchi et al. [2018] Ishibuchi et al. [2018] Ishibuchi et al. [2018] Ishibuchi et al. [2018] Balandat et al. [2020] Balandat et al. [2020] Balandat et al. [2020] Balandat et al. [2020] Balandat et al. [2020]

For TS-TCH [Paria et al., 2020], similarly to ParEGO, we employ random scalarisations by using the augmented Tchebycheff scalarisation function. After converting to single-objective optimisation problem, Thompson sampling is used as the acquisition function. We draw a sample from the joint posterior over a discrete set of 1000d points sampled from a scrambled Sobol sequence, suggested by Daulton et al. [2020]. In the batch setting, q distinct weight vectors are sampled, and the acquisition function is optimised sequentially for each. For EHVI [Daulton et al., 2020] and its noisy variant NEHVI [Daulton et al., 2021], the reference point is predefined, following the practice in Daulton et al. [2020], Yang et al. [2019a]. Table 4 lists the reference points used for each problem in the experimental evaluation. The logarithmic variants of both acquisition functions are employed, as recommended in Ament et al. [2023]. In the batch setting, the sequential greedy optimisation strategy is adopted. It is worth noting that in this study, EHVI and NEHVI are evaluated only on problems with 3 and 5 objectives, as the acquisition optimisation wall time becomes prohibitively high when the number of objectives increases to 10 (see Table 23). For C-EHVI [Gaudrie et al., 2018, 2020], the Kalai-Smorodinsky equilibrium is used to determine the reference point of HV. The logarithmic variant of the acquisition functions is employed, as recommended in Ament et al. [2023]. For joint entropy search (JES), we use S = 10 Monte Carlo samples and p = 10 number of Pareto optimal points, according to Tu et al. [2022]. In the batch setting, the sequential greedy optimisation strategy is adopted. For all the problems, we normalise the input variables and standardise the objective values before training Gaussian processes. We assume an independent surrogate model for each objective, using a constant mean function and a Matérn 5/2 ARD kernel. We optimise all acquisition functions by using the L-BFGS-B, with up to 200 iterations. 6

TS and BoTorch multi-objective tutorial EHVI, NEHVI, and BoTorch multi-objective tutorial 8 Botorch JES implementation 7

arXiv Template

E.3

A P REPRINT

Problem Details

The details of the benchmark problems and real-world problems are given in the following. DTLZ1.

DTLZ1 is a scalable benchmark problem from Deb et al. [2005], which is defined as:

1 x1 x2 · · · xm−1 (1 + g(xm )), 2 1 f2 (x) = x1 x2 · · · (1 − xm−1 )(1 + g(xm )), 2 .. . 1 fm−1 (x) = x1 (1 − x2 )(1 + g(xm )), 2 1 fm (x) = (1 − x1 )(1 + g(xm )), 2 s.t. 0 ≤ xi ≤ 1, for i = 1, 2, . . . , d. P where g(xm ) = 100 |xm | + xi ∈xm (xi − 0.5)2 − cos(20π(xi − 0.5)) and xm is the last d − m + 1 variables. f1 (x) =

DTLZ2.

DTLZ2 is a scalable benchmark problem from Deb et al. [2005], which is defined as: f1 (x) = (1 + g(xm )) cos(x1 π/2) · · · cos(xm−2 π/2) cos(xm−1 π/2), f2 (x) = (1 + g(xm )) cos(x1 π/2) · · · cos(xm−2 π/2) sin(xm−1 π/2), f3 (x) = (1 + g(xm )) cos(x1 π/2) · · · sin(xm−2 π/2), .. . fm (x) = (1 + g(xm )) sin(x1 π/2), s.t. 0 ≤ xi ≤ 1, for i = 1, 2, . . . , d.

where g(xm ) =

Inverted DTLZ1.

xi ∈xm (xi − 0.5)

and xm is the last d − m + 1 variables.

Inverted DTLZ1 is a variant of DTLZ1 [Jain and Deb, 2013a], which is defined as: fi (x) = 0.5 · (1 + g(xm )) − fiDTLZ1 (x), i = 1, . . . , m

where g(xm ) is the same function as used in DTLZ1, fiDTLZ1 (x) denotes the ith objective of the original DTLZ1 formulation. Inverted DTLZ2.

Inverted DTLZ2 is a variant of DTLZ2 [Jain and Deb, 2013b], which is defined as: fi (x) = 1 + g(xm ) − fiDTLZ2 (x), i = 1, . . . , m

where g(xm ) is the same function as used in DTLZ2, fiDTLZ2 (x) denotes the ith objective of the original DTLZ2 formulation. Convex DTLZ2.

Convex DTLZ2 is a variant of DTLZ2 [Deb and Jain, 2013], which is defined as: fi (x) = (fiDTLZ2 (x))4 , i = 1, . . . , m − 1 DTLZ2 fm (x) = (fm (x))2

where fiDTLZ2 (x) denotes the ith objective of the original DTLZ2 formulation. This problem convert the original concave problem to convex problem. Scaled DTLZ2.

Scaled DTLZ2 is a variant of DTLZ2, which is defined as: fi (x) = 2i−1 · fiDTLZ2 (x), i = 1, . . . , m

where fiDTLZ2 (x) denotes the ith objective of the original DTLZ2 formulation. This benchmark problem is used to see whether an algorithm can deal with problems with different scales of different objectives. 23

arXiv Template

A P REPRINT

DTLZ3–DTLZ7. DTLZ3–DTLZ7 are scalable multi-objective benchmark problems. Their mathematical formulations are provided in Deb et al. [2005]. According to the original paper [Deb et al., 2005], the dimensionality d of DTLZ1 and its variant (i.e., inverted DTLZ1) is m + 4, the dimensionality d of DTLZ2-6 and their variants (i.e., inverted DTLZ2, convex DTLZ2, and scaled DTLZ2) is m + 9, and the dimensionality d of DTLZ7 is m + 19. Car Side Impact Design. The car side-impact problem aims to minimise vehicle weight while satisfying safety constraints related to occupant injury and structural response [Jain and Deb, 2013a]. It involves m = 4 objectives with d = 7 variables, which are based on a surrogate model that is fit to data collected from a simulator. The mathematical formulations are given as follows:

f1 (x) = 1.98 + 4.9x1 + 6.67x2 + 6.98x3 + 4.01x4 + 1.78x5 + 10−5 x6 + 2.73x7 f2 (x) = 4.72 − 0.5x4 − 0.19x2 x3 f3 (x) = 0.5 (VMBP (x) + VFD (x)) f4 (x) = −

10 X

max (gi (x), 0)

i=1

where the constraint functions gi (x) are defined as: g1 (x) = 1 − 1.16 + 0.3717x2 x4 + 0.0092928x3 g2 (x) = 0.32 − 0.261 + 0.0159x1 x2 + 0.06486x1 + 0.019x2 x7 − 0.0144x3 x5 − 0.0154464x6 g3 (x) = 0.32 − 0.214 − 0.00817x5 + 0.045195x1 + 0.0135168x1 − 0.03099x2 x6 + 0.018x2 x7 − 0.007176x3 − 0.023223x3 + 0.00364x5 x6 + 0.018x22

g4 (x) = 0.32 − 0.74 + 0.61x2 + 0.031296x3 + 0.031872x7 − 0.227x22 g5 (x) = 32 − 28.98 − 3.818x3 + 4.2x1 x2 − 1.27296x6 + 2.68065x7 g6 (x) = 32 − 33.86 − 2.95x3 + 5.057x1 x2 + 3.795x2 + 3.4431x7 − 1.45728 g7 (x) = 32 − 46.36 + 9.9x2 + 4.4505x1 g8 (x) = 4 − f2 (x) g9 (x) = 9.9 − VMBP (x) g10 (x) = 15.7 − VFD (x) with volume terms defined as: VMBP (x) = 10.58 − 0.674x1 x2 − 0.67275x2 VFD (x) = 16.45 − 0.489x3 x7 − 0.8435x6 x7 . The search space is defined as: x1 ∈ [0.5, 1.5], x2 ∈ [0.45, 1.35], x3 , x4 ∈ [0.5, 1.5], x5 ∈ [0.875, 2.625], x6 , x7 ∈ [0.4, 1.2].

Car Cab Design. This vehicle performance optimisation problem involves m = 9 objectives with d = 7 variables, relating to aspects such as car roominess, fuel economy, acceleration time, and road noise at various speeds [Deb and Jain, 2013]. The problem includes 7 decision variables and 4 stochastic variables, which are based on a surrogate model that is fit to data collected from a simulator, defined as: 24

arXiv Template

A P REPRINT

f1 (x) = 1.98 + 4.9x1 + 6.67x2 + 6.98x3 + 4.01x4 + 1.75x5 + 10−5 x6 + 2.73x7 f2 (x) = [1.16 − 0.3717x2 x4 − 0.00931x2 x10 − 0.484x3 x9 + 0.01343x6 x10 ]+ 0.261 − 0.0159x1 x2 − 0.188x1 x8 − 0.019x2 x7 + 0.0144x3 x5 + 0.8757x5 x10 1 f3 (x) = 0.32 + 0.08045x6 x9 + 0.00139x8 x11 + 0.00001575x10 x11 +    0.214 + 0.00817x5 − 0.131x1 x8 − 0.0704x1 x9 + 0.03099x2 x6 − 0.018x2 x7  1  + 0.0208x x + 0.121x x − 0.00364x x + 0.0007715x x  f4 (x) =  3 8 3 9 5 6 5 10   0.32 − 0.0005354x6 x10 + 0.00121x8 x11 + 0.00184x9 x10 − 0.018x22 + 0.74 − 0.61x2 − 0.163x3 x8 + 0.001232x3 x10 − 0.166x7 x9 + 0.227x22 f5 (x) = 0.32 +    28.98 + 3.818x3 − 4.2x1 x2 + 0.0207x5 x10 + 6.63x6 x9 − 7.77x7 x8 + 0.32x9 x10 1 1 f6 (x) =  ·  + 33.86 + 2.95x3 + 0.1792x10 − 5.057x1 x2 − 11x2 x8 − 0.0215x5 x10 − 9.98x7 x8  32 3 + 22x8 x9 + 46.36 − 9.9x2 − 12.9x1 x8 + 0.1107x3 x10 + 2 4.72 − 0.5x4 − 0.19x2 x3 − 0.0122x4 x10 + 0.009325x6 x10 + 0.000191x11 f7 (x) = 4.0 + 10.58 − 0.674x1 x2 − 1.95x2 x8 + 0.02054x3 x10 − 0.0198x4 x10 + 0.028x6 x10 f8 (x) = 9.9 + 2 16.45 − 0.489x3 x7 − 0.843x5 x6 + 0.0432x9 x10 − 0.0556x9 x11 − 0.000786x11 f9 (x) = 15.7 + where [·]+ denotes max(0, ·) and the search space is defined as: x1 ∈ [0.5, 1.5], x2 ∈ [0.45, 1.35], x3 , x4 ∈ [0.5, 1.5], x5 ∈ [0.875, 2.625], x6 , x7 ∈ [0.4, 1.2]. The four stochastic variables are defined as: x8 ∼ N (0.345, 0.0062 ),

x9 ∼ N (0.192, 0.0062 ),

x10 , x11 ∼ N (0, 102 ).

arXiv Template

Additional Experimental Results

F.1

Noiseless Cases

A P REPRINT

In this section, we present the results on the six noiseless problems (i.e., DTLZ1 and DTLZ2 along with their four variants) with 3 and 10 objectives. Tables 5, 6 and 7 show the distance-based metric (log distance), the HV of the best solution (in terms of its HV value) and the HV of all evaluated solutions obtained by the SPMO and the peer methods on the six noiseless problems, respectively. Figures 6, 7 and 8 present the violin plots, illustrating the distributions of the corresponding results reported in Tables 5, 6, and 7, respectively. In addition, Figure 9 presents the trajectories of the distance metric obtained by each method on the noiseless problems with 3 and 10 objectives. We also give the results of the problems DTLZ3–DLTZ7 with 5 objectives, shown in Tables 8, 9 and 10. Table 5: Results of the distance-based metric (log distance) obtained by the SPMO and the peer methods on the noiseless problems with 3 objectives (top) and 10 objectives (bottom) on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 (3) Mean (Std)

DTLZ2 (3) Mean (Std)

Inverted DTLZ1 (3) Mean (Std)

Inverted DTLZ2 (3) Mean (Std)

Convex DTLZ2 (3) Mean (Std)

Scaled DTLZ2 (3) Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI C-EHVI JES SPMO Method

3.8e+0 (3.3e–1)+ 3.7e+0 (2.0e–1)+ 3.9e+0 (2.6e–1)+ 3.6e+0 (1.4e–1)+ 3.6e+0 (1.2e–1)+ 3.5e+0 (1.1e–1)+ 3.3e+0 (2.9e–1) DTLZ1 (10) Mean (Std)

2.3e–1 (5.1e–2)+ 4.7e–3 (3.1e–3)+ 1.4e–1 (7.1e–2)+ 4.6e–3 (2.5e–3)+ 1.6e–3 (1.7e–3)+ 5.5e–3 (4.2e–3)+ 2.2e–4 (1.2e–4) DTLZ2 (10) Mean (Std)

4.4e+0 (4.0e–1)+ 3.5e+0 (5.7e–1)+ 4.3e+0 (3.6e–1)+ 4.0e+0 (2.1e–1)+ 4.0e+0 (2.6e–1)+ 4.1e+0 (1.7e–1)+ 2.8e+0 (6.6e–1) Inverted DTLZ1 (10) Mean (Std)

5.7e–2 (5.3e–2)+ -3.0e–1 (7.7e–3)+ -1.1e–1 (2.4e–2)+ -3.0e–1 (4.3e–3)+ -2.8e–1 (2.9e–2)+ -3.0e–1 (4.8e–3)+ -3.1e–1 (1.3e–5) Inverted DTLZ2 (10) Mean (Std)

-7.7e–2 (1.7e–1)+ -9.3e–1 (1.2e–1)+ -3.4e–1 (1.3e–1)+ -9.1e–1 (1.1e–1)+ -4.5e–1 (2.4e–1)+ -9.1e–1 (9.0e–2)+ -1.2e+0 (1.4e–2) Convex DTLZ2 (10) Mean (Std)

2.2e–1 (4.8e–2)+ 5.2e–3 (4.7e–3)+ 1.5e–1 (3.6e–2)+ 1.8e–2 (6.0e–2)+ 3.2e–3 (5.5e–3)+ 9.6e–3 (8.7e–3)+ 3.8e–5 (1.6e–5) Scaled DTLZ2 (10) Mean (Std)

6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0

Sobol ParEGO TS-TCH C-EHVI JES SPMO

3.8e+0 (2.4e–1)+ 3.7e+0 (2.6e–1)+ 3.8e+0 (3.9e–1)+ 3.7e+0 (1.5e–1)+ 3.4e+0 (2.1e–1)+ 2.8e+0 (5.5e–1)

2.4e–1 (4.1e–2)+ 1.2e–1 (8.3e–2)+ 2.1e–1 (2.3e–2)+ 6.4e–3 (5.1e–3)+ 1.5e–1 (6.6e–2)+ 1.3e–3 (2.6e–3)

5.2e+0 (2.5e–1)+ 3.5e+0 (4.7e–1)∼ 5.3e+0 (3.7e–1)+ 3.9e+0 (5.1e–1)+ 4.5e+0 (5.3e–1)+ 3.4e+0 (5.1e–1)

1.2e+0 (5.4e–2)+ 8.6e–1 (1.9e–2)+ 1.0e+0 (1.4e–2)+ 8.4e–1 (1.6e–2)+ 8.6e–1 (2.0e–2)+ 7.8e–1 (2.4e–2)

-4.5e–1 (2.2e–1)+ -1.8e+0 (2.3e–1)+ -6.7e–1 (2.5e–1)+ -5.5e–1 (3.2e–1)+ -1.7e+0 (3.0e–1)+ -3.2e+0 (1.5e–1)

2.4e–1 (5.1e–2)+ 1.1e–1 (6.3e–2)+ 2.1e–1 (3.6e–2)+ 5.5e–2 (7.5e–2)+ 1.2e–1 (8.5e–2)+ 1.7e–2 (3.7e–2)

Sum up +/∼/− 6/ 0/ 0 5/ 1/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0

Table 6: The HV of the best solution (in terms of its HV value) obtained by the proposed SPMO and the peer methods on the noiseless problems with 3 objectives (top) and 10 objectives (bottom) on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 (3) Mean (Std)

DTLZ2 (3) Mean (Std)

Inverted DTLZ1 (3) Mean (Std)

Inverted DTLZ2 (3) Mean (Std)

Convex DTLZ2 (3) Mean (Std)

Scaled DTLZ2 (3) Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI C-EHVI JES SPMO Method

5.3e+7 (3.3e+6)+ 5.6e+7 (1.5e+6)+ 5.3e+7 (3.1e+6)+ 5.7e+7 (7.7e+5)+ 5.7e+7 (8.3e+5)+ 5.6e+7 (8.2e+5)+ 5.8e+7 (1.7e+6) DTLZ1 (10) Mean (Std)

3.3e–2 (2.1e–2)+ 1.6e–1 (7.0e–3)+ 5.0e–2 (3.1e–2)+ 1.6e–1 (4.6e–3)+ 1.3e–1 (1.3e–2)+ 1.5e–1 (8.3e–3)+ 1.7e–1 (5.3e–3) DTLZ2 (10) Mean (Std)

4.4e+7 (5.7e+6)+ 5.5e+7 (4.2e+6)+ 4.6e+7 (5.2e+6)+ 5.1e+7 (2.1e+6)+ 5.1e+7 (2.4e+6)+ 5.0e+7 (2.0e+6)+ 5.9e+7 (3.5e+6) Inverted DTLZ1 (10) Mean (Std)

1.1e–1 (2.5e–2)+ 3.1e–1 (3.7e–3)+ 2.0e–1 (1.2e–2)+ 3.1e–1 (2.2e–3)+ 3.0e–1 (1.4e–2)+ 3.0e–1 (2.4e–3)+ 3.1e–1 (6.4e–6) Inverted DTLZ2 (10) Mean (Std)

1.9e–1 (1.1e–1)+ 7.2e–1 (5.1e–2)+ 3.7e–1 (8.4e–2)+ 6.8e–1 (4.6e–2)+ 4.5e–1 (1.4e–1)+ 7.1e–1 (3.6e–2)+ 7.9e–1 (5.1e–3) Convex DTLZ2 (10) Mean (Std)

3.7e–2 (1.8e–2)+ 1.6e–1 (9.5e–3)− 4.8e–2 (2.5e–2)+ 1.3e–1 (3.4e–2)− 1.4e–1 (1.7e–2)− 1.5e–1 (8.1e–3)− 1.2e–1 (7.5e–3) Scaled DTLZ2 (10) Mean (Std)

6/ 0/ 0 5/ 0/ 1 6/ 0/ 0 5/ 0/ 1 5/ 0/ 1 5/ 0/ 1

Sobol ParEGO TS-TCH C-EHVI JES SPMO

8.7e+25 (3.8e+24)+ 9.2e+25 (3.3e+24)+ 8.8e+25 (5.1e+24)+ 9.2e+25 (1.0e+24)+ 9.3e+25 (1.9e+24)+ 9.7e+25 (3.6e+24)

4.7e–2 (2.5e–2)+ 1.2e–1 (9.6e–2)+ 4.3e–2 (4.0e–2)+ 2.6e–1 (3.2e–2)+ 9.9e–2 (7.2e–2)+ 2.9e–1 (3.4e–2)

2.4e+25 (9.4e+24)+ 7.8e+25 (1.1e+25)∼ 2.2e+25 (1.3e+25)+ 7.0e+25 (1.3e+25)+ 4.9e+25 (1.8e+25)+ 8.2e+25 (9.3e+24)

5.9e-11 (3.2e-10)+ 2.2e–5 (1.4e–5)+ 1.6e–8 (2.6e–8)+ 3.4e–5 (1.4e–5)+ 2.1e–5 (1.1e–5)+ 1.3e–4 (5.1e–5)

7.7e–1 (2.4e–1)+ 1.9e+0 (1.0e–1)+ 1.0e+0 (2.5e–1)+ 9.1e–1 (3.3e–1)+ 1.8e+0 (1.3e–1)+ 2.3e+0 (3.5e–2)

5.0e–2 (2.1e–2)+ 1.3e–1 (8.3e–2)+ 4.1e–2 (2.5e–2)+ 2.0e–1 (9.0e–2)+ 1.3e–1 (9.7e–2)+ 2.6e–1 (4.6e–2)

Sum up +/∼/− 6/ 0/ 0 5/ 1/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0

arXiv Template

A P REPRINT

Table 7: The HV of all the solutions obtained by the proposed SPMO and the peer methods on the noiseless problems with 3 objectives (top) and 10 objectives (bottom) on 30 independent runs, respectively. The method with the best mean is highlighted in bold. The symbols “+”, “∼”, and “−” indicate that a method is statistically worse than, equivalent to, and better than SPMO, respectively. Method

DTLZ1 (3) Mean (Std)

DTLZ2 (3) Mean (Std)

Inverted DTLZ1 (3) Mean (Std)

Inverted DTLZ2 (3) Mean (Std)

Convex DTLZ2 (3) Mean (Std)

Scaled DTLZ2 (3) Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI C-EHVI JES SPMO Method

6.3e+7 (2.8e+5)∼ 6.3e+7 (1.1e+6)− 6.2e+7 (7.8e+5)∼ 6.4e+7 (1.0e+5)− 6.3e+7 (6.8e+5)− 6.4e+7 (5.0e+4)− 6.2e+7 (1.0e+6) DTLZ1 (10) Mean (Std)

4.9e–2 (2.7e–2)+ 5.4e–1 (5.2e–2)− 7.5e–2 (4.2e–2)+ 6.4e–1 (2.2e–2)− 3.3e–1 (5.4e–2)+ 5.6e–1 (5.3e–2)− 4.6e–1 (6.3e–2) DTLZ2 (10) Mean (Std)

5.2e+7 (2.9e+6)+ 5.8e+7 (2.6e+6)+ 5.2e+7 (3.3e+6)+ 5.8e+7 (2.4e+6)+ 5.7e+7 (3.2e+6)+ 6.1e+7 (1.6e+6)∼ 6.1e+7 (1.5e+6) Inverted DTLZ1 (10) Mean (Std)

2.1e–1 (2.6e–2)+ 6.7e–1 (1.1e–2)− 4.1e–1 (1.2e–2)− 7.0e–1 (3.9e–3)− 3.2e–1 (3.5e–2)+ 6.6e–1 (1.1e–2)− 3.7e–1 (1.5e–2) Inverted DTLZ2 (10) Mean (Std)

2.5e–1 (1.3e–1)+ 1.0e+0 (9.0e–2)∼ 5.4e–1 (1.1e–1)+ 1.0e+0 (6.9e–2)∼ 5.7e–1 (2.0e–1)+ 1.0e+0 (6.3e–2)∼ 1.0e+0 (4.2e–2) Convex DTLZ2 (10) Mean (Std)

6.2e–2 (2.5e–2)+ 5.4e–1 (7.5e–2)− 8.8e–2 (4.6e–2)+ 2.6e–1 (7.7e–2)− 3.3e–1 (6.5e–2)− 5.0e–1 (8.5e–2)− 1.6e–1 (5.2e–2) Scaled DTLZ2 (10) Mean (Std)

5/ 1/ 0 1/ 1/ 4 4/ 1/ 1 1/ 1/ 4 4/ 0/ 2 0/ 2/ 4

Sobol ParEGO TS-TCH C-EHVI JES SPMO

2.4e+28 (1.3e+28)− 9.0e+26 (9.8e+26)+ 2.9e+28 (2.5e+28)− 1.0e+26 (1.6e+24)+ 8.8e+26 (1.3e+27)+ 1.2e+28 (1.5e+28)

3.6e–1 (1.7e–1)+ 1.7e+0 (3.5e+0)+ 1.8e–1 (2.0e–1)+ 5.2e–1 (1.2e–1)+ 6.3e–1 (1.1e+0)+ 5.2e+1 (3.3e+1)

5.3e+25 (5.9e+25)+ 3.6e+27 (1.3e+28)+ 4.7e+25 (4.4e+25)+ 7.1e+25 (1.3e+25)+ 5.7e+27 (1.9e+28)+ 2.3e+28 (5.3e+28)

5.9e-11 (3.2e-10)+ 7.2e–4 (4.9e–4)+ 1.9e–8 (3.3e–8)+ 6.6e–5 (2.6e–5)+ 7.5e–4 (3.9e–4)+ 1.4e–1 (1.5e–1)

5.2e+0 (3.2e+0)+ 2.7e+2 (2.2e+2)+ 2.8e+1 (2.4e+1)+ 1.3e+0 (4.0e–1)+ 2.8e+2 (1.9e+2)+ 6.3e+3 (8.0e+3)

3.6e–1 (1.6e–1)+ 6.8e–1 (7.6e–1)+ 1.5e–1 (1.0e–1)+ 3.8e–1 (2.2e–1)+ 1.5e+0 (2.2e+0)+ 5.2e+1 (6.5e+1)

DTLZ2 (3 obj) 1e 1 (5) DTLZ1

DTLZ1 (3 obj)

Log distance

2 0 2

4.4

Inverted DTLZ1 (3 obj) 4

4.0

3.8

3.6 DTLZ2 (3 obj) 1e 1 Inverted 0.0

3.2

0.5

3 2 1

Convex DTLZ2 (3 obj)

3.4

1e 1 Scaled DTLZ2 (3 obj) 3 2 1

75 1.0 100 125 150 175 200

Number of evaluations

Sobol

ParEGO

5/ 0/ 1 6/ 0/ 0 5/ 0/ 1 6/ 0/ 0 6/ 0/ 0

4.2

Log distance

4.5 4.0 3.5 3.0 2.5

Sum up +/∼/−

TS-TCH

EHVI

C-EHVI

JES

SPMO (ours)

(a) Violin plots of the distance-based metric (log distance) obtained by each method on problems with 3 objectives. 5

DTLZ1 (10 obj)

Log distance

4.2

Log distance

4.4

Log distance

4.0 Inverted DTLZ2 (10 obj) 3.8 1.2 1.0

1e 1 3

0 1 0 1

3.4

Sobol

DTLZ1 (3)

3.6

0.8

DTLZ2 (10 obj)

Convex DTLZ2 (10 obj)

503

100

150

ParEGO

TS-TCH

C-EHVI

Number of evaluations

6 5 4 3 2

Inverted DTLZ1 (10 obj)

1e 1 Scaled DTLZ2 (10 obj) 4 3 2 1 200 0 1 JES SPMO (ours)

(b) Violin plots of the distance-based metric (log distance) obtained by each method on problems with 10 objectives.

Figure 6: Violin plots of the distance-based metric (log distance) obtained by the proposed SPMO and the peer methods on the noiseless problems with 3 and 10 objectives. Each violin represents the distribution of the distance-based metric obtained by a method over 30 independent runs.

arXiv Template

DTLZ1 (3 obj)

Log HV (single point)

1e1

4.4

1.78

4.2

1.77

4.0

1.76

Log HV (single point)

DTLZ1 (5)DTLZ2 (3 obj) 4

1.76

1.74 1.72

Inverted 3.6DTLZ2 (3 obj)

Convex DTLZ2 (3 obj)

3.4

Scaled DTLZ2 (3 obj) 2

3.2

2.0

1e1 Inverted DTLZ1 (3 obj)

1.78

3.8

1.5

1.80

Log distance

1.79

A P REPRINT

2.5

25 Sobol

3 4

100 125 150 175 200

Number of evaluations

ParEGO

TS-TCH

EHVI

C-EHVI

JES

DTLZ1 (3) (10 obj) DTLZ2

SPMO (ours)

5.99 5.98 5.97

Log HV (single point)

5.96

DTLZ1 (10 obj)

4.4

6.0

5.9

4.2

5.8

4.0

Log distance

Log HV (single point)

(a) Violin plots of the HV of the best solution obtained by each method on problems with 3 objectives.

1e1

3.8

0.75 1e1 Inverted DTLZ2 (10 obj) 1.00 3.6 1.25 1.50 3.4 1.75 50 2.00

Sobol

1.0 0.5 0.0 0.5 1.0 1.5

5.7

Convex DTLZ2 (10 obj)

Scaled DTLZ2 (10 obj)

2 4

100

150

TS-TCH

C-EHVI

Number of evaluations

ParEGO

1e1 Inverted DTLZ1 (10 obj)

2006 8

JES

SPMO (ours)

(b) Violin plots of the HV of the best solution obtained by each method on problems with 10 objectives.

Figure 7: Violin plots of the HV of the best solution (in terms of its HV value) obtained by the proposed SPMO and the peer methods on the noiseless problems with 3 and 10 objectives. Each violin represents the distribution of HV values obtained by a method over 30 independent runs.

arXiv Template

DTLZ1 (5)DTLZ2 (3 obj) 1e 1

DTLZ1 (3 obj) 4.4

1e7 6.4

Log HV

Log distance

4.0

6.0

3.8

3.6

1e 1 Inverted DTLZ2 (3 obj)

Convex DTLZ2 (3 obj)

3.4

Log HV

6 4

0.5

Sobol

100 125 150 175 200

0.0 Number of evaluations

ParEGO

1e7 Inverted DTLZ1 (3 obj)

1e 1 Scaled DTLZ2 (3 obj)

1.0

3.2

6.5 6.0 5.5 5.0 4.5 4.0

4.2

6.2

A P REPRINT

TS-TCH

EHVI

C-EHVI

2 0

JES

DTLZ1 (3)(10 obj) DTLZ2

SPMO (ours)

(a) Violin plots of the HV of all evaluated solutions obtained by each method on problems with 3 objectives.

Log HV

6.6 6.4 6.2 6.0

DTLZ1 (10 obj)

4.4 4.2

Log distance

6.8

1e1

4.0 3.8

Log HV

0.0 0.5 1.0

5.0 2.5 0.0 2.5 5.0 7.5

1e1 Inverted DTLZ2 (10 obj)

3.6 3.4

1.5 2.0

Sobol

7.00 6.75 6.50 6.25 6.00 5.75

1e1 Convex DTLZ2 (10 obj) 1.00 0.75 0.50 0.25 50 0.00

Scaled DTLZ2 (10 obj) 5 0

100

150

TS-TCH

C-EHVI

Number of evaluations

ParEGO

1e1 Inverted DTLZ1 (10 obj)

2005 JES

SPMO (ours)

(b) Violin plots of the HV of all evaluated solutions obtained by each method on problems with 10 objectives.

Figure 8: Violin plots of the HV of all evaluated solutions obtained by the proposed SPMO and the peer methods on the noiseless problems with 3 objectives (top) and 10 objectives (bottom), respectively. Each violin represents the distribution of HV values obtained by a method over 30 independent runs.

arXiv Template

4.0

4.2

0.2

3.5

4.0

0.1

Log distance

0.3

Log distance

4.4

50 3.8 100

0.2

Log distance

DTLZ1 (5)DTLZ2 (3 obj)

DTLZ1 (3 obj)

4.5

150

200

0.0

Number of evaluations Inverted 3.6 DTLZ2 (3 obj)

0.0 0.2

A P REPRINT

100

5.0 4.5 4.0 3.5 3.0

150

200

Number of evaluations Convex DTLZ2 (3 obj)

3.4

0.0

3.2

0.5

Inverted DTLZ1 (3 obj)

100

150

200

Number of evaluations 0.3

Scaled DTLZ2 (3 obj)

0.2 0.1

1.0

150 150 175 200 200 0.0 50 100 150 200 50 100 25 15050 20075 100 50125 100 Number of evaluations of evaluations Number of evaluations Number of Number evaluations Sobol

ParEGO

TS-TCH

EHVI

C-EHVI

JES

SPMO (ours)

(a) Trajectories of the distance metric on the problems with 3 objectives.

DTLZ1 (10 obj)

Log distance

4.0

3.0 50

4.2

0.1

100 4.0 150

Log distance

1.0 0.8

100

Inverted DTLZ1 (10 obj)

0.2

Number of evaluations Inverted DTLZ2 3.8 (10 obj) 1.2

DTLZ2 DTLZ1 (3)(10 obj)

4.4

Log distance

3.5

0.3

200

0.0

100

150

Number of evaluations

200

Convex DTLZ2 (10 obj)

0.3

3.6

0.2

3.4

0.1

150

20050

Sobol

ParEGO

Number of evaluations

100

150

TS-TCH

C-EHVI

50 100 150 200 Number of evaluations Number of evaluations

JES

100

150

Number of evaluations

200

Scaled DTLZ2 (10 obj)

200 50 100 150 200 Number of evaluations SPMO (ours)

(b) Trajectories of the distance metric on the problems with 10 objectives.

Figure 9: Trajectories of the distance metric obtained by the SPMO and the peer methods on the noiseless problems with 3 and 10 objectives. Each coloured line represents the mean distance of the closest solution to the utopian point on 30 independent runs (after the initial Sobol samples, represented by the dashed grey line).

arXiv Template

A P REPRINT

Table 8: Results of the distance-based metric (log distance) obtained by obtained by the SPMO and the peer methods on DTLZ3–DTLZ7 with 5 objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ3 Mean (Std)

DTLZ4 Mean (Std)

DTLZ5 Mean (Std)

DTLZ6 Mean (Std)

DTLZ7 Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI JES SPMO

6.3e+0 (1.6e–1)+ 5.5e+0 (1.1e–1)∼ 6.3e+0 (1.6e–1)+ 5.4e+0 (4.6e–2)− 5.5e+0 (5.0e–2)∼ 5.5e+0 (1.8e–1)

2.7e–1 (6.0e–2)+ 2.3e–1 (7.1e–2)+ 2.6e–1 (5.0e–2)+ 8.7e–2 (9.4e–2)+ 2.3e–1 (8.1e–2)+ 1.2e–2 (1.9e–2)

2.5e–1 (4.1e–2)+ 4.0e–2 (3.1e–2)+ 2.3e–1 (5.9e–2)+ 3.0e–1 (6.0e–2)+ 5.2e–2 (5.8e–2)+ 3.2e–3 (5.0e–3)

2.2e+0 (2.1e–2)+ 4.3e–2 (1.6e–1)+ 2.2e+0 (1.6e–2)+ -9.4e–7 (8.4e–7)+ 2.1e–2 (1.1e–1)+ -1.9e–6 (5.6e–7)

3.1e+0 (5.0e–2)+ 2.2e+0 (8.0e–2)+ 3.1e+0 (4.8e–2)+ 1.7e+0 (8.5e–2)+ 2.2e+0 (5.8e–2)+ 1.7e+0 (1.1e–1)

5/ 0/ 0 4/ 1/ 0 5/ 0/ 0 4/ 0/ 1 4/ 1/ 0

Table 9: The HV of the best solution (in terms of its HV value) obtained by SPMO and the peer methods on the DTLZ3–DTLZ7 with 5 objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ3 Mean (Std)

DTLZ4 Mean (Std)

DTLZ5 Mean (Std)

DTLZ6 Mean (Std)

DTLZ7 Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI JES SPMO

9.2e+19 (1.1e+18)+ 9.7e+19 (5.2e+17)∼ 9.3e+19 (1.0e+18)+ 9.8e+19 (1.1e+17)− 9.8e+19 (1.4e+17)− 9.7e+19 (1.0e+18)

2.2e–2 (1.4e–2)+ 2.5e–2 (3.0e–2)+ 1.3e–2 (1.1e–2)+ 1.3e–1 (5.6e–2)+ 3.6e–2 (3.5e–2)+ 1.6e–1 (2.5e–2)

8.4e+4 (9.5e+2)∼ 8.9e+4 (1.4e+3)− 8.6e+4 (1.1e+3)− 8.3e+4 (1.6e+3)+ 8.9e+4 (1.5e+3)− 8.5e+4 (2.6e+3)

8.7e+3 (8.2e+2)+ 8.8e+4 (3.8e+3)− 9.3e+3 (9.9e+2)+ 9.0e+4 (1.3e–2)− 8.8e+4 (3.5e+3)− 8.3e+4 (2.9e+3)

-0.0e+0 (0.0e+0)+ 2.8e+5 (2.7e+4)+ -0.0e+0 (0.0e+0)+ 3.8e+5 (2.2e+4)∼ 2.8e+5 (2.0e+4)+ 3.8e+5 (2.8e+4)

4/ 1/ 0 2/ 1/ 2 4/ 0/ 1 2/ 1/ 2 2/ 0/ 3

Table 10: The HV of all evaluated solutions obtained by SPMO and the peer methods on DTLZ3–DTLZ7 with 5 objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ3 Mean (Std)

DTLZ4 Mean (Std)

DTLZ5 Mean (Std)

DTLZ6 Mean (Std)

DTLZ7 Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI JES SPMO

3.9e+20 (1.7e+20)∼ 2.1e+20 (2.1e+20)∼ 5.0e+20 (3.0e+20)− 1.7e+20 (1.5e+20)∼ 2.7e+20 (2.7e+20)∼ 4.2e+20 (7.0e+20)

4.9e–2 (2.4e–2)+ 5.3e–2 (7.4e–2)+ 2.1e–2 (1.8e–2)+ 9.3e–1 (7.1e–1)∼ 8.8e–2 (1.3e–1)+ 1.1e+0 (9.3e–1)

2.1e+5 (1.1e+5)∼ 2.4e+5 (8.7e+4)∼ 1.8e+5 (8.2e+4)+ 3.5e+5 (1.8e+5)− 3.6e+5 (3.5e+5)∼ 2.5e+5 (1.4e+5)

1.4e+5 (1.7e+4)+ 3.9e+5 (1.6e+5)− 1.8e+5 (4.2e+4)+ 4.2e+5 (1.2e+5)− 4.0e+5 (1.6e+5)− 3.1e+5 (1.9e+5)

-0.0e+0 (0.0e+0)+ 9.1e+5 (4.1e+5)+ -0.0e+0 (0.0e+0)+ 1.9e+6 (5.3e+5)∼ 1.0e+6 (4.3e+5)+ 2.3e+6 (1.3e+6)

3/ 2/ 0 2/ 2/ 1 4/ 0/ 1 0/ 3/ 2 2/ 2/ 1

arXiv Template

F.2

A P REPRINT

Noisy Cases

In this section, we present the results on the noisy problems. Tables 11, 12 and 13 show the distance-based metric (log distance), the HV of the best solution (in terms of its HV value) and the HV of all evaluated solutions obtained by the SPMO and the peer methods, respectively. Figures 10, 11 and 12 present the violin plots, illustrating the distributions of the corresponding results reported in Tables 11, 12 and 13, respectively. Table 11: Results of the distance-based metric (log distance) obtained by the SPMO and the peer methods on the noisy problems with 5 objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

Sobol NParEGO TS-TCH NEHVI JES SPMO

3.7e+0 (4.2e–1)+ 3.6e+0 (3.1e–1)+ 3.8e+0 (3.5e–1)+ 3.5e+0 (1.4e–1)+ 3.4e+0 (1.2e–1)+ 3.0e+0 (6.0e–1)

1.9e–1 (5.6e–2)+ 1.0e–1 (8.2e–2)+ 2.0e–1 (4.2e–2)+ -1.1e–1 (8.1e–2)+ 1.3e–1 (1.2e–1)+ -1.8e–1 (5.8e–2)

4.7e+0 (3.4e–1)+ 3.1e+0 (3.8e–1)+ 4.9e+0 (3.2e–1)+ 4.0e+0 (5.8e–1)+ 4.5e+0 (1.1e–1)+ 2.9e+0 (4.8e–1)

6.1e–1 (6.1e–2)+ 1.8e–1 (5.0e–2)+ 4.6e–1 (5.5e–2)+ 1.5e–1 (3.3e–2)+ 1.9e–1 (6.5e–2)+ 5.8e–2 (6.8e–2)

-2.9e–1 (2.2e–1)+ -1.3e+0 (2.7e–1)+ -4.7e–1 (2.3e–1)+ -1.3e+0 (3.0e–1)+ -6.9e–1 (3.5e–1)+ -2.1e+0 (3.3e–1)

2.2e–1 (4.4e–2)+ 5.2e–2 (5.9e–2)+ 2.1e–1 (5.0e–2)+ 2.9e–1 (9.5e–2)+ 1.0e–1 (8.7e–2)+ -1.9e–1 (1.8e–1)

6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0

Table 12: The HV of the best solution (in terms of its HV value) obtained by SPMO and the peer methods on the noisy problems with 5 objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

Sobol NParEGO TS-TCH NEHVI JES SPMO

8.6e+12 (5.2e+11)+ 9.1e+12 (2.5e+11)+ 8.5e+12 (5.1e+11)+ 9.1e+12 (1.5e+11)+ 9.0e+12 (1.3e+11)+ 9.4e+12 (4.4e+11)

5.1e–2 (3.0e–2)+ 1.0e–1 (4.9e–2)+ 5.0e–2 (1.4e–2)+ 2.8e–1 (7.6e–2)+ 8.7e–2 (7.2e–2)+ 3.2e–1 (4.9e–2)

5.4e+12 (1.2e+12)+ 9.0e+12 (4.1e+11)+ 4.8e+12 (1.2e+12)+ 7.4e+12 (1.2e+12)+ 6.1e+12 (3.2e+11)+ 9.2e+12 (4.2e+11)

8.2e–4 (1.5e–3)+ 5.8e–2 (1.6e–2)+ 8.1e–3 (5.4e–3)+ 6.7e–2 (1.0e–2)+ 5.8e–2 (2.0e–2)+ 1.0e–1 (2.7e–2)

4.0e–1 (1.8e–1)+ 1.2e+0 (2.0e–1)+ 5.3e–1 (1.7e–1)+ 1.2e+0 (2.1e–1)+ 7.8e–1 (3.1e–1)+ 1.8e+0 (2.2e–1)

3.8e–2 (1.9e–2)+ 1.3e–1 (4.0e–2)+ 3.1e–2 (1.6e–2)+ 1.2e–2 (1.0e–2)+ 1.0e–1 (5.7e–2)+ 3.8e–1 (1.6e–1)

6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0

Table 13: The HV of all the solutions obtained by the proposed SPMO and the peer methods on the noisy problems with five objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼”, and “−” indicate that a method is statistically worse than, equivalent to, and better than SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

Sobol NParEGO TS-TCH NEHVI JES SPMO

4.5e+13 (1.4e+13)+ 5.8e+13 (2.3e+13)∼ 5.1e+13 (1.3e+13)∼ 7.2e+13 (2.1e+13)∼ 7.5e+13 (2.3e+13)∼ 6.2e+13 (2.9e+13)

2.1e–1 (8.4e–2)+ 4.3e–1 (3.9e–1)+ 1.9e–1 (8.3e–2)+ 3.1e+0 (8.9e–1)− 4.3e–1 (5.0e–1)+ 2.0e+0 (5.6e–1)

8.2e+12 (4.7e+12)+ 1.7e+13 (9.6e+12)+ 6.9e+12 (2.6e+12)+ 1.5e+13 (1.7e+13)+ 3.5e+13 (2.0e+13)∼ 4.4e+13 (3.6e+13)

8.5e–4 (1.5e–3)+ 5.1e–1 (1.3e–1)+ 3.2e–2 (1.8e–2)+ 7.7e–1 (2.0e–1)+ 5.7e–1 (1.5e–1)+ 9.6e–1 (3.1e–1)

9.3e–1 (5.1e–1)+ 1.1e+1 (4.8e+0)+ 2.2e+0 (9.8e–1)+ 9.5e+0 (3.8e+0)+ 4.6e+0 (2.8e+0)+ 1.4e+1 (4.4e+0)

1.6e–1 (5.7e–2)+ 7.7e–1 (6.5e–1)+ 9.0e–2 (6.5e–2)+ 1.5e–2 (1.3e–2)+ 4.3e–1 (3.3e–1)+ 2.9e+0 (1.8e+0)

6/ 0/ 0 5/ 1/ 0 5/ 1/ 0 4/ 1/ 1 4/ 2/ 0

arXiv Template DTLZ1 (5 obj)

Log distance

DTLZ2 (5 obj)

1e 1

Inverted DTLZ1 (5 obj) 5

5.0

Convex DTLZ2 (5 obj)

Log distance

Inverted DTLZ2 (5 obj) 8 1e 1 4.5

6 4 2 0 2

Inverted DTLZ1 (5)

5.5

Log distance

A P REPRINT

5.0

4.0

2.5

3.5

0.0

3.0 25

Sobol

1e 1 Scaled DTLZ2 (5 obj)

2.5 5.0

100 3 125 150 175 200

NParEGO

TS-TCH

NEHVI

JES

SPMO (ours)

Figure 10: Violin plots of the distance-based metric (log distance) obtained by the six methods on the noisy problems with five objectives. Each violin represents the distribution of the distance-based metric obtained by a method over 30 independent runs. Log HV (single point)

3.00

1e1

DTLZ1 (5 obj)

2.99

3.00

2.98 2.96

5.0

0.50 0.75 1.00

Log distance

1e1 Inverted DTLZ2 (5 obj) 0.25

4.5

4.0

3.5

2.85

Convex DTLZ2 (5 obj)

Scaled DTLZ2 (5 obj) 0 2 4 6 8

3.0

1.25

2.90

Inverted DTLZ1 (5)

5.5

Sobol

100 125 150 175 200

NParEGO

1e1 Inverted DTLZ1 (5 obj)

2.95

2.97

Log HV (single point)

DTLZ2 (5 obj)

TS-TCH

NEHVI

JES

SPMO (ours)

Figure 11: Violin plots of the HV of the best solution (in terms of its HV value) obtained by the six methods on the noisy problems with five objectives. Each violin represents the distribution of maximum HV values obtained by a method over 30 independent runs. 1e14

DTLZ1 (5 obj)

Log HV

0.5

Inverted DTLZ1 (5)

5.5

5.0

Log HV

1.5 1.0 0.5 0.0

Inverted DTLZ2 (5 obj)

Log distance

2.0

2.0 1.5 1.0 0.5 0.0

1.0 0.0

DTLZ2 (5 obj)

1.5

4.5

4.0

3.5

3.0 25

Sobol

1e1 Convex DTLZ2 (5 obj)

8 6 4 2 0 2

100 0 125 150 175 200

NParEGO

TS-TCH

NEHVI

JES

1e14 Inverted DTLZ1 (5 obj)

Scaled DTLZ2 (5 obj)

SPMO (ours)

Figure 12: Violin plots of the HV of all evaluated solutions obtained by the six methods on the noisy problems with five objectives. Each violin represents the distribution of maximum HV values obtained by a method over 30 independent runs.

arXiv Template

F.3

A P REPRINT

Batch Setting

We compare the proposed SPMO with the peer methods in the batch setting where the batch size q is set to 5 (a commonly used value [Lin et al., 2022]). Tables 14, 15 and 16 show the distance-based metric (log distance), the HV of the best solution (in terms of its HV value) and the HV of all evaluated solutions obtained by the SPMO and the peer methods, respectively. Figures 13, 14 and 15 present the violin plots, illustrating the distributions of the corresponding results reported in Tables 14, 15 and 16, respectively. Table 14: Results of the distance-based metric (log distance) obtained by the SPMO and the five peer methods with a batch size q = 5 on the problems with 5 objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI JES SPMO

3.7e+0 (3.0e–1)+ 3.4e+0 (1.9e–1)+ 3.9e+0 (2.2e–1)+ 3.5e+0 (9.6e–2)+ 3.4e+0 (1.3e–1)+ 3.1e+0 (3.0e–1)

2.4e–1 (4.6e–2)+ 6.7e–2 (8.1e–2)+ 2.0e–1 (4.1e–2)+ 9.3e–3 (3.5e–3)+ 1.1e–1 (8.0e–2)+ 9.0e–4 (8.3e–4)

4.8e+0 (3.1e–1)+ 3.2e+0 (5.4e–1)+ 4.8e+0 (3.2e–1)+ 4.0e+0 (4.4e–1)+ 4.5e+0 (1.7e–1)+ 2.9e+0 (4.9e–1)

6.0e–1 (5.7e–2)+ 2.5e–1 (1.3e–2)+ 4.6e–1 (1.2e–2)+ 2.3e–1 (5.5e–3)+ 2.6e–1 (2.0e–2)+ 2.1e–1 (5.0e–5)

-3.1e–1 (2.5e–1)+ -1.7e+0 (2.0e–1)+ -7.1e–1 (2.1e–1)+ -1.4e+0 (2.1e–1)+ -1.0e+0 (4.5e–1)+ -2.1e+0 (7.7e–3)

2.3e–1 (5.1e–2)+ 1.3e–1 (8.9e–2)+ 2.6e–1 (3.8e–2)+ 3.1e–1 (6.0e–2)+ 8.8e–2 (9.3e–2)+ 1.7e–4 (1.2e–4)

6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 6/ 0/ 0

Table 15: The HV of the best solution (in terms of its HV value) obtained by SPMO and the five peer methods with a batch size q = 5 on the problems with 5 objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI JES SPMO

8.7e+12 (3.7e+11)+ 9.1e+12 (1.8e+11)+ 8.4e+12 (3.8e+11)+ 9.0e+12 (9.5e+10)+ 9.0e+12 (1.2e+11)+ 9.3e+12 (2.8e+11)

3.5e–2 (1.5e–2)+ 1.3e–1 (5.6e–2)+ 3.2e–2 (1.3e–2)+ 1.9e–1 (7.2e–3)∼ 1.0e–1 (3.9e–2)+ 1.9e–1 (1.4e–2)

5.1e+12 (1.1e+12)+ 8.8e+12 (7.0e+11)+ 4.9e+12 (1.2e+12)+ 7.5e+12 (9.7e+11)+ 6.2e+12 (5.0e+11)+ 9.2e+12 (4.8e+11)

8.7e–4 (1.5e–3)+ 4.0e–2 (2.9e–3)+ 7.1e–3 (1.1e–3)+ 4.5e–2 (1.4e–3)+ 3.9e–2 (4.5e–3)+ 4.9e–2 (1.2e–5)

3.8e–1 (1.7e–1)+ 1.1e+0 (6.2e–2)+ 6.5e–1 (1.4e–1)+ 1.0e+0 (9.7e–2)+ 8.5e–1 (2.4e–1)+ 1.3e+0 (5.0e–3)

3.8e–2 (2.1e–2)+ 8.7e–2 (5.3e–2)+ 2.8e–2 (1.5e–2)+ 1.5e–2 (1.4e–2)+ 1.1e–1 (5.5e–2)+ 1.6e–1 (1.6e–2)

6/ 0/ 0 6/ 0/ 0 6/ 0/ 0 5/ 1/ 0 6/ 0/ 0

Table 16: The HV of all the solutions obtained by the six methods on the problems with five objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼”, and “−” indicate that a method is statistically worse than, equivalent to, and better than SPMO, respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

Sobol ParEGO TS-TCH EHVI JES SPMO

4.4e+13 (1.4e+13)∼ 2.9e+13 (1.2e+13)∼ 4.5e+13 (1.3e+13)∼ 2.9e+13 (1.8e+13)∼ 3.1e+13 (2.8e+13)∼ 3.7e+13 (2.2e+13)

1.4e–1 (5.6e–2)+ 6.5e–1 (5.3e–1)+ 8.9e–2 (5.3e–2)+ 2.1e+0 (5.9e–1)− 6.9e–1 (5.2e–1)+ 1.5e+0 (7.0e–1)

7.7e+12 (2.9e+12)+ 1.5e+13 (7.8e+12)+ 6.8e+12 (1.8e+12)+ 1.9e+13 (1.7e+13)+ 7.4e+13 (9.1e+13)− 3.7e+13 (3.5e+13)

9.5e–4 (1.6e–3)+ 5.0e–1 (9.6e–2)∼ 5.0e–2 (1.7e–2)+ 6.1e–1 (1.5e–1)∼ 4.7e–1 (1.1e–1)∼ 5.9e–1 (3.3e–1)

9.5e–1 (4.6e–1)+ 9.6e+0 (4.8e+0)+ 3.2e+0 (1.7e+0)+ 7.4e+0 (2.2e+0)+ 5.7e+0 (4.6e+0)+ 2.6e+1 (1.7e+1)

1.6e–1 (6.6e–2)+ 4.2e–1 (3.9e–1)+ 7.4e–2 (5.2e–2)+ 2.0e–2 (1.7e–2)+ 7.5e–1 (6.3e–1)+ 1.9e+0 (1.6e+0)

5/ 1/ 0 4/ 2/ 0 5/ 1/ 0 3/ 2/ 1 3/ 2/ 1

arXiv Template DTLZ2(3) (5 obj) 1e 1 DTLZ1

DTLZ1 (5 obj)

4.4

Log distance

4 3

Log distance

4.2

4.0

Log distance

3.6

3.4

4 2

Sobol

Inverted DTLZ1 (5 obj)

3.8(5 obj) 1e 1 Inverted DTLZ2

A P REPRINT

Convex DTLZ2 (5 obj) 0.0 0.5 1.0 1.5 50 2.0

100

150

TS-TCH

EHVI

Number of evaluations

ParEGO

1e 1 Scaled DTLZ2 (5 obj) 4 3 2 1 0 200 1

JES

SPMO (ours)

Figure 13: Violin plots of the distance-based metric (log distance) obtained by the proposed SPMO and the peer methods on the problems with five objectives. Each violin represents the distribution of the distance-based metric obtained by a method over 30 independent runs. Log HV (single point)

3.00 1e1 2.99

DTLZ2 (5 obj)

4.4

4.2

4.0

Log distance

2.98

DTLZ1 (3)

DTLZ1 (5 obj)

2.97 2.96

3.8

Log HV (single point)

1e1 Inverted DTLZ2 (5 obj) 0.4 0.6 0.8 1.0 1.2

3.6 3.4

2.90 2.85

100

150

TS-TCH

EHVI

Number of evaluations

ParEGO

1e1 Inverted DTLZ1 (5 obj)

2.95

Convex DTLZ2 (5 obj) 0

Sobol

3.00

1e1 Scaled DTLZ2 (5 obj) 0.0 0.2 0.4 0.6 0.8 200 1.0

JES

SPMO (ours)

Figure 14: Violin plots of the HV of the best solution (in terms of its HV value) obtained by the six methods with a batch size q = 5 on the problems with five objectives. Each violin represents the distribution of maximum HV values obtained by a method over 30 independent runs.

Log HV

1.0 0.5 0.0

Log HV

1.25 1.00 0.75 0.50 0.25 0.00

1e14

DTLZ1 (3)

DTLZ1 (5 obj)

DTLZ2 (5 obj)

4.4

4.2

Log distance

1.5

1.00 0.75 0.50 0.25 0.00

4.0

3.8

Inverted DTLZ2 (5 obj)

1e1 Convex DTLZ2 (5 obj)

Scaled DTLZ2 (5 obj)

3.6

3.4

Sobol

1e14 Inverted DTLZ1 (5 obj)

100

150

TS-TCH

EHVI

0 Number of evaluations

ParEGO

200 0 JES

SPMO (ours)

Figure 15: Violin plots of the HV of all evaluated solutions obtained by the six methods on the problems with five objectives. Each violin represents the distribution of maximum HV values obtained by a method over 30 independent runs.

arXiv Template

F.4

A P REPRINT

Sensitivity analysis

In this section, we conduct a sensitivity analysis to assess the effect of different utopian points. We consider three different settings. The first one is slightly better than the ideal point, i.e., with a difference of 0.01, the second is fairly better than the ideal point (i.e. 0.1), and the last one is significantly better than the ideal point (i.e. 1.0). Tables 17, 18 and 19 show the distance-based metric (log distance), the HV of the best solution (in terms of its HV value), and the HV of all evaluated solutions obtained by the SPMO with four different utopian points, respectively. Figures 16, 17 and 18 present the violin plots, illustrating the distributions of the corresponding results reported in Tables 17, 18 and 19, respectively. Table 17: Results of the distance-based metric (log distance) obtained by the SPMO with four different utopian points on the problems with five objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than SPMO (current), respectively. Method SPMO_0.01 SPMO_0.1 SPMO_1.0 SPMO

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

3.1e+0 (4.9e–1)∼ 2.9e+0 (3.6e–1)− 2.8e+0 (6.1e–1)∼ 3.1e+0 (3.0e–1)

6.2e–4 (1.4e–3)− 3.0e–4 (2.3e–4)− 3.3e–4 (2.0e–4)− 9.0e–4 (8.3e–4)

2.8e+0 (4.2e–1)∼ 3.0e+0 (5.3e–1)∼ 2.9e+0 (5.0e–1)∼ 2.9e+0 (4.9e–1)

2.1e–1 (3.6e–5)∼ 2.1e–1 (4.6e–5)∼ 2.1e–1 (3.4e–5)∼ 2.1e–1 (5.0e–5)

-2.1e+0 (1.6e–2)∼ -2.1e+0 (2.6e–2)∼ -2.1e+0 (1.6e–2)∼ -2.1e+0 (7.7e–3)

6.4e–5 (3.9e–5)− 5.3e–5 (2.4e–5)− 5.7e–5 (2.4e–5)− 1.7e–4 (1.2e–4)

0/ 4/ 2 0/ 3/ 3 0/ 4/ 2

Table 18: The HV of the best solution (in terms of its HV value) obtained by the proposed SPMO with four different utopian points on the problems with five objectives on 30 independent runs. The method exhibiting the best mean is highlighted in bold. The symbols “+”, “∼”, and “−” denote that a method is statistically worse than, equivalent to, or better than SPMO (current), respectively. Method SPMO_0.01 SPMO_0.1 SPMO_1.0 SPMO

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

9.4e+12 (3.7e+11)∼ 9.5e+12 (2.7e+11)− 9.5e+12 (3.8e+11)− 9.3e+12 (2.8e+11)

1.9e–1 (1.0e–2)∼ 1.9e–1 (1.6e–2)∼ 1.9e–1 (1.3e–2)∼ 1.9e–1 (1.4e–2)

9.3e+12 (3.8e+11)∼ 9.0e+12 (5.6e+11)∼ 9.1e+12 (4.2e+11)∼ 9.2e+12 (4.8e+11)

4.9e–2 (9.0e–6)∼ 4.9e–2 (1.1e–5)∼ 4.9e–2 (8.5e–6)∼ 4.9e–2 (1.2e–5)

1.3e+0 (7.8e–3)∼ 1.2e+0 (1.1e–2)∼ 1.3e+0 (7.8e–3)∼ 1.3e+0 (5.0e–3)

1.6e–1 (1.4e–2)∼ 1.6e–1 (1.0e–2)∼ 1.6e–1 (1.4e–2)∼ 1.6e–1 (1.6e–2)

0/ 6/ 0 0/ 5/ 1 0/ 5/ 1

Table 19: The HV of all the solutions obtained by the SPMO with four different utopian points on the problems with five objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼”, and “−” indicate that a method is statistically worse than, equivalent to, and better than SPMO (current), respectively. Method SPMO_0.01 SPMO_0.1 SPMO_1.0 SPMO

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

3.1e+13 (1.7e+13)∼ 3.5e+13 (1.5e+13)∼ 3.3e+13 (1.4e+13)∼ 3.7e+13 (2.2e+13)

2.7e+0 (1.2e+0)− 2.7e+0 (1.2e+0)− 3.0e+0 (1.6e+0)− 1.5e+0 (7.0e–1)

3.5e+13 (3.7e+13)∼ 2.6e+13 (1.7e+13)∼ 3.2e+13 (1.9e+13)∼ 3.7e+13 (3.5e+13)

6.5e–1 (7.0e–1)∼ 6.1e–1 (4.5e–1)∼ 7.1e–1 (6.3e–1)∼ 5.9e–1 (3.3e–1)

2.5e+1 (1.8e+1)∼ 2.7e+1 (1.7e+1)∼ 2.9e+1 (2.2e+1)∼ 2.6e+1 (1.7e+1)

2.2e+0 (1.5e+0)∼ 1.8e+0 (1.0e+0)∼ 2.1e+0 (1.5e+0)∼ 1.9e+0 (1.6e+0)

0/ 5/ 1 0/ 5/ 1 0/ 5/ 1

arXiv Template

Inverted DTLZ1 (5)

5.5

DTLZ1 (5 obj)

4.5

2 1

4.0

4 3 2 1

Inverted DTLZ2 (5 obj)

1e 1

Log distance

2.122 2.121

3.5

Convex DTLZ2 (5 obj)

1.95

8 4

2.05

2.119

SPMO_0.01

1e 4 Scaled DTLZ2 (5 obj)

2.00

3.0

2.120

Inverted DTLZ1 (5 obj)

8 6 4 2 0

Log distance

DTLZ2 (5 obj)

1e 3

5.0

A P REPRINT

50 2.10 75

100 125 150 1750 200

Number of evaluations SPMO_0.1

SPMO_1.0

SPMO (current)

Figure 16: Violin plots of the distance-based metric (log distance) obtained by the four methods on the problems with five objectives. Each violin represents the distribution of the distance-based metric obtained by a method over 30 independent runs.

Inverted DTLZ1 (5)

5.5 2.99 2.98

Log HV (single point)

2.97

3.0145 3.0150

DTLZ1 (5 obj)

4.5 4.0 3.5

2.4

3.0160

SPMO_0.01

1e1 Inverted DTLZ1 (5 obj)

2.99 2.98 2.97 2.96

Scaled DTLZ2 (5 obj)

1e 1

2.2

3.0 25

3.00

Convex DTLZ2 (5 obj)

Inverted DTLZ2 (5 obj)

3.0155

DTLZ2 (5 obj)

1.5 1.6 1.7 1.8 1.9 2.0

5.0

Log distance

Log HV (single point)

3.00 1e1

2.0

100 125 150

1.6 1.7 1.8 1.9 1752.0

1.8 Number of evaluations

SPMO_0.1

SPMO_1.0

200

SPMO (current)

Figure 17: Violin plots of the HV of the best solution (in terms of its HV value) obtained by the four methods on the problems with five objectives. Each violin represents the distribution of maximum HV values obtained by a method over 30 independent runs.

Inverted DTLZ1 (5)

5.5

Log HV

DTLZ1 (5 obj)

1.00

5.0 4.5

0.50

4.0

Log HV

DTLZ2 (5 obj)

1e14 Inverted DTLZ1 (5 obj) 2.0 1.5 1.0 0.5 0.0

0.00

0.25

Inverted DTLZ2 (5 obj) 4 3 2 1 0

1e1

0.75

Log distance

1e14 1.00 0.75 0.50 0.25 0.00

1e2 Convex DTLZ2 (5 obj)

3.5

1.0

3.0

0.5

25 SPMO_0.01

50 0.0 75

100 125 150

Scaled DTLZ2 (5 obj) 8 6 4 2 1750

Number of evaluations SPMO_0.1

SPMO_1.0

200

SPMO (current)

Figure 18: Violin plots of the HV of all evaluated solutions obtained by the four methods the problems with five objectives. Each violin represents the distribution of maximum HV values obtained by a method over 30 independent runs.

arXiv Template

F.5

A P REPRINT

Comparison of Single-Point Metrics within SPMO

In the proposed SPMO framework, we employ a distance metric (i.e., the distance of a solution to the utopian point). However, different metrics can be adopted provided that it can reflect the quality of a solution in achieving a good trade1 1 off between objectives, such as the weighted sum and Tchebycheff scalarisation with the same weights ( m ,..., m ), where m denotes the number of objectives. Here, we compare these three variants of SPMO, i.e., SPMOdist , SPMOT ch , and SPMOws . Tables 20, 21 and 22 show the distance-based metric (log distance), the HV of the best solution (in terms of its HV value) and the HV of all evaluated solutions obtained by the SPMO with three different single-point metrics, respectively. Figures 19, 20 and 21 present the violin plots, illustrating the distributions of the corresponding results reported in Tables 20, 21 and 22, respectively. Table 20: Results of the distance-based metric (log distance) obtained by the SPMO using three different single-point metrics on the problems with five objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼” and “−” indicate that the method is statistically worse than, equivalent to and better than our SPMO (i.e., SPMOdist ), respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

SPMOT ch SPMOws SPMO

3.5e+0 (3.2e–1)+ 3.6e+0 (1.3e–1)+ 3.1e+0 (3.0e–1)

3.5e–5 (4.3e–5)− 2.1e–4 (1.5e–4)− 9.0e–4 (8.3e–4)

2.9e+0 (6.5e–1)∼ 3.1e+0 (5.4e–1)∼ 2.9e+0 (4.9e–1)

3.1e–1 (3.3e–2)+ 2.2e–1 (2.7e–3)+ 2.1e–1 (5.0e–5)

-1.8e+0 (1.7e–1)+ -1.2e+0 (2.7e–1)+ -2.1e+0 (7.7e–3)

2.9e–5 (1.9e–5)− 3.2e–4 (3.4e–4)+ 1.7e–4 (1.2e–4)

3/ 1/ 2 4/ 1/ 1

Table 21: The HV of the best solution (in terms of its HV value) obtained by the proposed SPMO using three different single-point metrics on the problems with five objectives on 30 independent runs. The method exhibiting the best mean is highlighted in bold. The symbols “+”, “∼”, and “−” denote that a method is statistically worse than, equivalent to, or better than SPMO (i.e., SPMOdist ), respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

SPMOT ch SPMOws SPMO

9.1e+12 (3.0e+11)+ 9.1e+12 (1.3e+11)+ 9.3e+12 (2.8e+11)

1.6e–1 (1.9e–2)+ 1.7e–1 (1.4e–2)+ 1.9e–1 (1.4e–2)

9.1e+12 (5.7e+11)∼ 9.0e+12 (6.5e+11)∼ 9.2e+12 (4.8e+11)

2.7e–2 (5.9e–3)+ 4.8e–2 (6.7e–4)+ 4.9e–2 (1.2e–5)

1.2e+0 (4.2e–2)+ 9.7e–1 (1.2e–1)+ 1.3e+0 (5.0e–3)

1.6e–1 (2.2e–2)∼ 1.7e–1 (1.5e–2)∼ 1.6e–1 (1.6e–2)

4/ 2/ 0 4/ 2/ 0

Table 22: The HV of all the solutions obtained by the SPMO using three different single-point metrics on the problems with five objectives on 30 independent runs. The method with the best mean is highlighted in bold. The symbols “+”, “∼”, and “−” indicate that a method is statistically worse than, equivalent to, and better than SPMO (i.e., SPMOdist ), respectively. Method

DTLZ1 Mean (Std)

DTLZ2 Mean (Std)

Inverted DTLZ1 Mean (Std)

Inverted DTLZ2 Mean (Std)

Convex DTLZ2 Mean (Std)

Scaled DTLZ2 Mean (Std)

Sum up +/∼/−

SPMOT ch SPMOws SPMO

1.5e+13 (6.2e+12)+ 1.3e+13 (3.1e+12)+ 3.7e+13 (2.2e+13)

4.9e+0 (7.6e+0)− 1.6e+0 (7.5e–1)∼ 1.5e+0 (7.0e–1)

2.1e+13 (1.5e+13)+ 2.4e+13 (1.8e+13)+ 3.7e+13 (3.5e+13)

2.1e+0 (2.6e+0)− 7.4e–1 (7.5e–1)∼ 5.9e–1 (3.3e–1)

1.3e+1 (7.4e+0)+ 6.4e+0 (3.8e+0)+ 2.6e+1 (1.7e+1)

3.3e+0 (3.4e+0)− 1.2e+0 (3.5e–1)∼ 1.9e+0 (1.6e+0)

3/ 0/ 3 3/ 3/ 0

arXiv Template DTLZ1 (5 obj) 4.0

1e 3

Log distance

3.0

Log distance

4.5 1e 1 Inverted DTLZ2 4.0 (5 obj)

Log distance

4.0 3.0

2.0 1.5 1.0 0.5 0.0

1.0 1.5

2.0

Convex DTLZ2 (5 obj)

3.0

2.5

0.5

3.5

4.5

2.0

Inverted DTLZ1 (5 obj)

5.0

2.5

DTLZ2 (5 obj)

Inverted DTLZ1 (5) 4

5.5

3.5

A P REPRINT

1e 3 Scaled DTLZ2 (5 obj)

100 125 150 175 200

SPMO_Tch

SPMO_ws

SPMO (SPMO_dist)

Figure 19: Violin plots of the distance-based metric (log distance) obtained by the SPMO using three different singlepoint metrics on the problems with five objectives. Each violin represents the distribution of the distance-based metric obtained by a method over 30 independent runs. 1e1

2.990 2.985 2.980

Log HV (single point)

2.975 3.0 3.5 4.0

DTLZ1 (5 obj)

DTLZ2 (5 obj)

Inverted DTLZ1 (5)

5.5

1.8

5.0

2.0

4.5

Inverted 4.0 DTLZ2 (5 obj)

1e 1 Convex DTLZ2 (5 obj) 2

3.5

Scaled DTLZ2 (5 obj) 1.8

2.0

1e1 Inverted DTLZ1 (5 obj)

1.6

3.0

4.5

3.00 2.99 2.98 2.97 2.96 2.95

1.6

Log distance

Log HV (single point)

2.995

756 100 125 150 175 2002.2

SPMO_Tch

SPMO_ws

SPMO (SPMO_dist)

Figure 20: Violin plots of the HV of the best solution (in terms of its HV value) obtained by the SPMO using three different single-point metrics on the problems with 5 objectives. Each violin represents the distribution of the singlepoint HV obtained by a method over 30 independent runs. 1e14

DTLZ1 (5 obj)

1e1

5.5

Log HV

0.5 0.0

1e14 Inverted DTLZ1 (5 obj) 2.0 1.5 1.0 0.5 0.0

5.0

2 0

4.5

DTLZ2 (5 obj) 1e1 Inverted 4.0 1.0

DTLZ2 (5 obj)

Inverted DTLZ1 (5)

Log distance

Log HV

1.00 0.75 0.50 0.25 0.00

3.5 3.0 25

SPMO_Tch

8 6 4 2 0

1e1 Convex DTLZ2 (5 obj)

100 125 150 175

SPMO_ws

1e1 Scaled DTLZ2 (5 obj) 2.0 1.5 1.0 0.5 2000.0

SPMO (SPMO_dist)

Figure 21: Violin plots of the HV of all evaluated solutions obtained by the SPMO using three different single-point metrics on the problems with five objectives. Each violin represents the distribution of maximum HV values obtained by a method over 30 independent runs.

arXiv Template

F.6

A P REPRINT

Acquisition Wall Time

Table 23 presents the mean acquisition optimisation wall time of the eight methods. As shown, when the number of objectives is 3 or 5, the time of all the methods is acceptable with a maximum of 98 seconds. As the number of objectives increases to 10, hypervolume-based methods (i.e., EHVI and NEHVI) become very expensive (taking about half an hour and more than 3 hours, respectively). The proposed SPMO method shows high computational efficiency, achieving the lowest time requirement in four out of the six instances. Table 23: Mean acquisition optimisation wall time in seconds based on the 2(d + 1) initial Sobol samples on DTLZ1 problems with m = 3, 5, 10 objectives, where d = m + 4, over 30 runs. Experiments are conducted using a CPU (Intel Xeon CPU Platinum 8360Y @ 2.40 GHz) and a GPU (NVIDIA A100). Note that N/A means the wall time of NEHVI on DTLZ1 with 10 objectives exceeds 3 hours.

Device\Method

ParEGO

NParEGO

TS-TCH

EHVI

NEHVI

C-EHVI

JES

SPMO (ours)

CPU (3 obj) GPU (3 obj) CPU (5 obj) GPU (5 obj) CPU (10 obj) GPU (10 obj)

3.46 5.19 2.59 9.02 8.76 41.33

3.49 3.94 2.14 10.33 5.97 30.30

12.29 14.74 28.32 22.29 87.44 64.78

4.23 3.05 36.03 10.32 1134.35 2426.04

6.15 5.73 97.98 63.90 N/A N/A

10.23 9.48 29.89 23.89 77.94 86.45

9.24 14.02 25.16 25.80 388.37 676.31

2.40 2.24 2.58 5.31 6.95 24.05

Applicability of the Proposed SPMO

In conventional multi-objective optimisation, algorithms are designed to approximate the entire Pareto front, so that a decision-maker can later select a preferred solution based on their own preferences. This is the ideal situation, as a well-represented Pareto front provides the most comprehensive view of the possible trade-offs [Jiang and Li, 2025b]. However, under tight evaluation budgets - especially when many objectives are involved - it is often unrealistic, if not impossible, to obtain a good approximation of the Pareto front. In such settings, the decision-maker may benefit more from receiving a well-balanced solution that is close to the Pareto front, rather than a well-distributed solution set far from the front. The proposed SPMO framework is designed for this purpose: instead of spreading search effort across the whole front, it directs the optimisation towards a well-balanced solution. With a focused search effort, the obtained solution is often closer to the front, and thus has a higher likelihood of being selected by the decision-maker. It is worth pointing out that if the optimisation problem under consideration is extremely simple (e.g., smooth, unimodal landscape with a very limited search space), on which finding a good representation of the entire Pareto front is possible under tight budgets, then our approach may not be desirable. Moreover, exploring the Pareto front can help decision-makers better understand the optimisation problem and facilitate the elicitation or refinement of their preferences. In such cases, our framework is not applicable.

Extensions

The preceding discussion has addressed multi-objective optimisation problems in which evaluations are performed either sequentially or in batches, under both noiseless and noisy scenarios. However, not all multi-objective settings conform to these scenarios. To accommodate a broader class of problems, we propose several extensions that enable the methodology to handle more optimisation scenarios. High-Dimensional Bayesian Optimisation (HDBO). High-dimensional black-box optimisation problems are highly challenging and frequently encountered in a wide range of applications. The dimensionality may range from tens to a billion [González-Duque et al., 2024, Papenmeier et al., 2023, Santoni et al., 2024, Wang et al., 2016, Hoang et al., 2025]. To tackle such optimisation problems, various HDBO methods have been proposed [Binois and Wycoff, 2022, Chen et al., 2024, Nayebi et al., 2019, Wang et al., 2018, 2016, Xu et al., 2025]. They can loosely be categorised into four classes [Santoni et al., 2024], i.e., variable selection [Eriksson and Jankowiak, 2021], additive models [Delbridge et al., 2020, Han et al., 2021, Wang et al., 2018, Ziomek and Ammar, 2023], embeddings [Antonov et al., 2022, Letham et al., 2020, Raponi et al., 2020], and trust regions [Daulton et al., 2022b, Diouane et al., 2023, Eriksson et al., 2019]. However, recent studies show that standard Gaussian processes without the above techniques can perform well in high-dimensional spaces [Hvarfner et al., 2024, Papenmeier et al., 2025, Xu et al., 2025] and suggest that the main issue in high-dimensional BO is the gradient vanishing. Our work can be naturally extended to the high-dimensional setting by mitigating the gradient vanishing issue [Papenmeier et al., 2025, Xu et al., 2025]. 40

arXiv Template

A P REPRINT

Multi-Fidelity Bayesian Optimisation (MFBO). In many real-world optimisation scenarios, the evaluation is often available at multiple fidelity levels, where increasing fidelity typically leads to improved accuracy at the expense of higher computational cost. Many MFBO methods have been proposed to tackle such optimisation problems [Belakaria et al., 2020b, Kandasamy et al., 2017, Li et al., 2020, Moss et al., 2021, Song et al., 2019, Takeno et al., 2020, Wu et al., 2020, Zhang et al., 2017]. Our proposed SPMO can be potentially extended to the multi-fidelity setting by integrating prior techniques, e.g., building multiple surrogate models of different levels of fidelity.

Record · ID 6006 · SHA-256 6835368de9b5f3e5
Conceptio Open Knowledge Archive — every document is proof-bundled with source, license, and retrieval metadata.