Scalarizing Functions in Bayesian Multiobjective Optimization

TL;DR

This paper reviews and compares 15 scalarizing functions within Bayesian multiobjective optimization, analyzing their impact on model performance and evaluation efficiency across various benchmark problems.

Contribution

It provides a comprehensive review and empirical comparison of scalarizing functions in Bayesian multiobjective optimization, highlighting their influence on optimization quality and efficiency.

Findings

Scalarizing functions significantly affect the number of evaluations needed.

Different functions perform variably depending on the problem and objectives.

Insights aid in selecting appropriate scalarizing functions for Bayesian optimization.

Abstract

Scalarizing functions have been widely used to convert a multiobjective optimization problem into a single objective optimization problem. However, their use in solving (computationally) expensive multi- and many-objective optimization problems in Bayesian multiobjective optimization is scarce. Scalarizing functions can play a crucial role on the quality and number of evaluations required when doing the optimization. In this article, we study and review 15 different scalarizing functions in the framework of Bayesian multiobjective optimization and build Gaussian process models (as surrogates, metamodels or emulators) on them. We use expected improvement as infill criterion (or acquisition function) to update the models. In particular, we compare different scalarizing functions and analyze their performance on several benchmark problems with different number of objectives to be optimized. The review and experiments on different functions provide useful insights when using and selecting a scalarizing function when using a Bayesian multiobjective optimization method.

Tables

Table 1. Table 1: Number of variables ( n 𝑛 n ) for WFG and DTLZ suites

	WFG			DTLZ
$k$	$d$	$l$	$n$	$n$
2	4	4	8	6
3	4	4	8	7
5	8	4	12	9
10	18	4	22	14

Table 2. Table 2: Parameters values used in different scalarizing functions

Scalarizing function	Parameter value
EWC	$p$ = 100
WPO	$p$ = 3
WN	$p$ = 0.5
ATCH	$\alpha$ = 0.0001
MTCH	$\alpha$ = 0.0001
PBI	$\theta$ = 5
IPBI	$\theta$ = 5
QPBI	$\theta$ = 1
APD	$\alpha$ = 2
HypI	$r$ (In table 3 for DTLZ problems)

Table 3. Table 3: Reference point in calculating hypevolumes in HypI and as performance measure:

Problem	reference point
DTLZ1	400 $\times$ (1,k)
DTLZ2	1.5 $\times$ (1,k)
DTLZ3	900 $\times$ (1,k)
DTLZ4	2 $\times$(1,k)
DTLZ5	1.5 $\times$ (1,k)
DTLZ6	6 $\times$ (1,k)
DTLZ7	5 $\times$ (1,k)

Problem	k	WS	EWC	WPO	WN	WPR	TCH	ATCH	MTCH
	2	42.26 (8.13)	46.28 (16.52)	37.95 (20.37)	41.54 (7.98)	38.38 (9.24)	41.15 (7.59)	38.86 (5.55)	38.98 (7.00)
DTLZ1	3	41.55 (5.51)	41.03 (14.44)	33.07 (7.59)	42.60 (5.81)	34.86 (5.83)	39.33 (3.66)	39.81 (2.91)	40.32 (2.86)
	5	38.40 (5.30)	34.56 (11.79)	29.82 (8.37)	38.85 (5.24)	37.24 (10.61)	36.92 (5.07)	38.88 (5.72)	38.23 (3.84)
	10	38.65 (12.18)	39.07 (15.04)	34.51 (9.55)	39.29 (13.78)	31.14 (11.38)	35.91 (12.20)	37.54 (10.52)	38.09 (9.53)
	2	0.19 (0.03)	0.14 (0.03)	0.05 (0.01)	0.19 (0.03)	0.19 (0.03)	0.04 (0.00)	0.04 (0.00)	(0.00)
DTLZ2	3	0.26 (0.02)	0.21 (0.02)	0.14 (0.01)	0.27 (0.02)	0.25 (0.01)	0.16 (0.01)	0.16 (0.01)	0.17 (0.02)
	5	0.41 (0.02)	0.39 (0.04)	0.32 (0.01)	0.40 (0.02)	0.37 (0.02)	0.40 (0.01)	0.40 (0.02)	0.40 (0.01)
	10	0.81 (0.04)	0.82 (0.03)	0.76 (0.02)	0.80 (0.05)	0.83 (0.04)	0.80 (0.03)	0.81 (0.04)	0.81 (0.03)
	2	93.74 (10.84)	123.41 (43.55)	94.69 (20.13)	88.45 (12.71)	(8.14)	88.54 (14.42)	89.96 (12.05)	92.26 (12.10)
DTLZ3	3	99.90 (8.77)	115.48 (20.65)	90.15 (25.95)	90.38 (8.50)	89.41 (14.54)	99.88 (6.58)	97.63 (8.99)	93.69 (12.77)
	5	93.10 (13.14)	137.71 (30.78)	105.19 (32.37)	100.31 (11.86)	93.92 (28.96)	96.68 (15.43)	95.11 (17.12)	95.96 (11.63)
	10	101.40 (15.91)	116.94 (30.53)	103.88 (22.65)	109.04 (25.48)	117.87 (32.94)	90.99 (21.19)	103.90 (14.91)	93.17 (11.40)
	2	0.34 (0.04)	0.41 (0.16)	0.34 (0.03)	0.38 (0.09)	0.45 (0.09)	0.32 (0.05)	0.33 (0.03)	0.33 (0.03)
DTLZ4	3	0.50 (0.05)	0.54 (0.08)	0.46 (0.03)	0.56 (0.05)	0.54 (0.08)	0.47 (0.03)	0.48 (0.04)	0.47 (0.04)
	5	0.64 (0.03)	0.67 (0.06)	0.57 (0.02)	0.64 (0.03)	0.69 (0.04)	0.64 (0.03)	0.65 (0.04)	0.65 (0.04)
	10	0.76 (0.03)	0.80 (0.03)	0.75 (0.03)	0.77 (0.02)	0.82 (0.02)	0.75 (0.03)	0.76 (0.02)	0.76 (0.03)
	2	0.21 (0.03)	0.14 (0.02)	0.05 (0.01)	0.21 (0.03)	0.20 (0.03)	0.04 (0.01)	0.04 (0.01)	0.05 (0.00)
DTLZ5	3	0.18 (0.02)	0.13 (0.02)	(0.01)	0.18 (0.02)	0.18 (0.02)	0.06 (0.01)	0.07 (0.01)	0.06 (0.01)
	5	0.18 (0.03)	0.12 (0.02)	0.04 (0.01)	0.18 (0.03)	0.17 (0.03)	0.09 (0.02)	0.08 (0.01)	0.09 (0.02)
	10	0.15 (0.02)	0.11 (0.01)	0.06 (0.01)	0.15 (0.02)	0.15 (0.02)	0.14 (0.02)	0.14 (0.03)	0.14 (0.03)
	2	0.41 (0.03)	3.67 (0.32)	(0.04)	0.40 (0.02)	0.62 (0.18)	0.36 (0.08)	0.37 (0.09)	0.39 (0.09)
DTLZ6	3	0.46 (0.07)	3.58 (0.24)	0.49 (0.14)	0.53 (0.11)	0.82 (0.19)	0.42 (0.10)	0.46 (0.07)	0.44 (0.06)
	5	0.74 (0.12)	3.71 (0.26)	1.27 (0.42)	0.75 (0.15)	1.19 (0.32)	0.85 (0.13)	0.75 (0.18)	0.88 (0.18)
	10	1.71 (0.50)	3.51 (0.33)	1.62 (0.46)	2.26 (0.59)	2.99 (0.50)	1.57 (0.39)	1.68 (0.29)	1.48 (0.32)
	2	0.22 (0.02)	3.27 (0.89)	0.30 (0.05)	0.24 (0.08)	0.35 (0.03)	0.25 (0.02)	0.25 (0.03)	0.25 (0.01)
DTLZ7	3	0.47 (0.05)	5.10 (1.16)	0.69 (0.13)	0.56 (0.04)	0.58 (0.04)	0.44 (0.03)	0.46 (0.03)	0.45 (0.03)
	5	0.98 (0.05)	8.16 (2.35)	2.67 (1.76)	0.98 (0.04)	0.99 (0.04)	0.92 (0.02)	0.92 (0.07)	0.94 (0.03)
	10	1.55 (0.08)	13.92 (3.39)	4.55 (3.09)	(0.05)	1.73 (0.14)	1.58 (0.07)	1.56 (0.07)	1.56 (0.06)

Table 5. Table 4: Mean IGD values and standard deviation (in parentheses) for DTLZ problems. The values statistically similar to the best one are in bold, the best value is encircled the worst value is underlined

Problem	k	PBI	IPBI	HypI	DomRank	MSD	QPBI	APD
	2	36.95 (6.20)	41.75 (12.41)	70.82 (24.42)	36.93 (10.38)	(12.12)	38.19 (11.76)	36.44 (4.26)
DTLZ1	3	34.63 (6.66)	43.24 (10.63)	60.13 (21.74)	37.79 (15.00)	32.37 (10.85)	(11.70)	34.58 (8.18)
	5	(8.23)	36.59 (12.28)	48.57 (16.87)	34.74 (10.43)	21.29 (13.64)	26.69 (5.78)	28.94 (10.56)
	10	31.60 (16.09)	36.73 (11.10)	39.49 (15.73)	39.97 (14.51)	34.94 (7.95)	35.80 (14.57)	(15.99)
	2	0.11 (0.01)	0.08 (0.02)	0.08 (0.02)	0.10 (0.03)	0.13 (0.05)	0.09 (0.02)	0.04 (0.01)
DTLZ2	3	0.15 (0.01)	0.26 (0.02)	(0.01)	0.18 (0.02)	0.16 (0.03)	0.15 (0.01)	0.14 (0.01)
	5	0.25 (0.01)	0.45 (0.01)	0.28 (0.02)	0.36 (0.01)	0.30 (0.06)	(0.00)	0.23 (0.00)
	10	0.83 (0.04)	0.85 (0.03)	(0.06)	0.77 (0.04)	0.75 (0.06)	0.81 (0.04)	0.74 (0.04)
	2	93.65 (14.09)	108.98 (20.84)	177.05 (64.60)	100.48 (36.21)	98.52 (13.61)	90.55 (11.43)	94.19 (13.95)
DTLZ3	3	91.46 (28.57)	98.01 (27.37)	169.47 (53.32)	112.14 (27.98)	78.94 (24.08)	(25.97)	104.18 (23.84)
	5	58.51 (39.42)	102.92 (24.46)	170.63 (55.24)	117.51 (39.97)	75.52 (28.84)	101.55 (11.68)	(46.26)
	10	(29.42)	108.02 (36.65)	126.97 (55.58)	118.07 (23.69)	102.90 (27.71)	98.89 (11.89)	99.95 (46.68)
	2	0.30 (0.06)	0.32 (0.05)	0.26 (0.08)	0.33 (0.05)	0.34 (0.03)	0.33 (0.07)	(0.08)
DTLZ4	3	0.45 (0.04)	0.49 (0.06)	0.50 (0.13)	0.50 (0.04)	0.41 (0.09)	0.42 (0.03)	(0.07)
	5	(0.02)	0.73 (0.09)	0.82 (0.09)	0.63 (0.05)	0.71 (0.08)	0.58 (0.04)	0.56 (0.02)
	10	(0.04)	0.88 (0.02)	0.91 (0.03)	0.78 (0.03)	0.76 (0.03)	0.72 (0.03)	0.78 (0.04)
	2	0.12 (0.01)	0.08 (0.01)	0.07 (0.01)	0.09 (0.02)	0.14 (0.06)	0.09 (0.01)	(0.01)
DTLZ5	3	0.05 (0.01)	0.14 (0.02)	0.05 (0.01)	0.09 (0.02)	0.08 (0.02)	0.04 (0.01)	0.06 (0.01)
	5	0.08 (0.02)	0.17 (0.02)	0.07 (0.01)	0.08 (0.02)	0.06 (0.02)	0.03 (0.01)	(0.00)
	10	0.13 (0.03)	0.15 (0.02)	0.10 (0.02)	0.09 (0.02)	0.08 (0.02)	0.09 (0.01)	(0.01)
	2	0.32 (0.10)	1.48 (0.24)	0.48 (0.22)	0.73 (0.51)	0.36 (0.08)	0.63 (0.28)	0.41 (0.18)
DTLZ6	3	(0.04)	1.53 (0.44)	1.50 (0.62)	2.69 (0.66)	1.68 (1.71)	0.57 (0.15)	0.24 (0.05)
	5	(0.04)	2.33 (0.53)	2.48 (0.93)	2.52 (0.76)	3.60 (0.51)	2.24 (0.38)	0.34 (0.05)
	10	(0.08)	3.55 (0.47)	3.20 (0.45)	3.26 (0.51)	3.61 (0.45)	2.43 (0.50)	1.01 (0.50)
	2	0.24 (0.02)	0.37 (0.00)	(0.02)	0.41 (0.11)	0.32 (0.21)	0.73 (0.13)	0.23 (0.09)
DTLZ7	3	0.81 (0.13)	0.85 (0.07)	(0.12)	0.73 (0.14)	0.63 (0.32)	1.25 (0.33)	0.67 (0.19)
	5	2.20 (0.36)	5.15 (2.82)	(0.13)	2.20 (0.62)	1.95 (0.93)	2.93 (0.23)	2.03 (0.40)
	10	3.84 (0.49)	13.56 (3.45)	6.86 (5.17)	7.89 (4.05)	11.42 (4.56)	8.28 (2.66)	6.03 (2.87)

Equations

$$\begin{array}{c}\text{mbox}minimize\left\{f_1(x),\dots ,f_k(x)\right\}\\ \text{mbox}subjecttox\in S\end{array}$$

$$g=\sum _{i=1}^k{w}_i{f}_i.$$

$$g=\sum _{i=1}^k\exp (p{w}_i-1)\exp (p{f}_i).$$

$$g=\sum _{i=1}^k{w}_i(f_i{)}^p.$$

$$g=(\sum _{i=1}^k{w}_i|f_i{|}^p{)}^{1/p}.$$

$$g=\prod _{i=1}^k(f_i{)}^{w_i}.$$

$$g=\underset{i}{\max }[w_i|f_i-z_i^{\ast }|],$$

$$g=\underset{i}{\max }[w_i|f_i-z_i^{\ast }|]+\alpha \sum _{i=1}^k|f_i-z_i^{\ast }|$$

$$g=\max_{i}\big{[}w_{i}(|f_{i}-z_{i}^{*}|+\alpha\sum_{i=1}^{k}|f_{i}-z_{i}^{*}|)\big{]}$$

$$g=d_1+\theta d_2,$$

$$g=\theta d_2-d_1,$$

$$g=d_1+\theta d_2\frac{d_2}{d^{\ast }},$$

$$d^{\ast }=\alpha \frac{1}{H}\frac{1}{k}\sum (z^{nadir}-z^{ideal})$$

$$g=1+P(\theta )\cdot \Vert \mathbf{f}\Vert ,$$

$$g=H(x\cup P_{k+1}).$$

$$g=1-\frac{rank(x)-1}{|X|-1}.$$

$$g=\min d(x^{\prime },x)$$

Full text

Keywords: metamodelling, machine learning, multiple criteria decision making, Pareto optimality, computational cost, Bayesian optimization

1 INTRODUCTION

Bayesian multiobjective optimization methods have been widely used for solving (computationally) expensive multiobjective optimization problems (MOPs) and in the last few years, many articles have been published. A survey of many Bayesian and non-Bayesian methods with their advantages and limitations is provided in [6, 29, 30, 1, 10, 19].

A generic framework of a Bayesian optimization method is shown in Figure 1, where the steps are mentioned inside squares.

In the first step, a set of samples is generated e.g. using a design of experiment technique (DOE) like  Latin hypercube sampling [41]. These samples are then evaluated with expensive objective functions in the second step. A termination criterion (typically maximum number of expensive function evaluations) is checked in the third step. If the termination criterion is met, nondominated solutions from all evaluated solutions (i.e. from expensive evaluations) are used as the final solutions. Otherwise, surrogate models are built using evaluated solutions in the fourth step. An optimizer e.g. an evolutionary algorithm is then used with the models in the fifth step to find promising samples by using an appropriate infilling criterion (or updating criterion or acquisition function). In the sixth step, a fixed number of samples generated with the optimizer is selected which are then evaluated with expensive objective functions.

All six steps mentioned above are important in the performance of a Bayesian optimization method. In this article, our focus is on the fourth step i.e. building or training surrogates. Once, the samples are evaluated with expensive objective functions in the second step, there can be different ways to build surrogates to be used in a Bayesian optimization method. The most common way is to build surrogate for each objective function [5, 47, 55, 37]. For instance, in [5, 55, 47], Gaussian process models are used for different objective functions. Another way is to build a surrogate for a scalarizing function after converting multiobjective optimization problem into a single objective optimization problem [33, 48], which is also the focus of this article. There are other approaches for building surrogates like classification of solutions into different ranks or classes [44, 38, 39, 46, 3]. For details on other ways of building surrogates, see [6, 1]. Note that scalarizing function can also be used after building models on each objective function e.g. as in [55, 34] to do a local search. In this article, we focus on building surrogate after converting a multiobjective optimization problem into a single objective by using a scalarizing function.

There are two main advantages to make a single surrogate in solving an expensive MOP. The first one is that only one surrogate is used in the solution process instead of multiple surrogates, which reduces the computational burden e.g. the training time especially in many-objective (usually more than three objectives) optimization problems. The second advantage is that one can use an infill criterion proposed for single-objective optimization problems which also reduces the computational complexity. In the literature, a little attention has been paid in using scalarizing functions for building surrogates on them and only a few studies exist in the literature. For instance, in [33], an augmented achievement scalarizing function (AASF) was used and in [48], hypervolume improvement, dominance rank and minimum signed distance were used (more details are provided in the Section 3).

It is worth important to be pointed out that several studies exist [45, 23, 16, 43, 42] on using scalarizing functions without using them in a Bayesian optimization method. In this article, we study 15 different scalarizing functions to build surrogates on them in solving (computationally) expensive MOPs and particularly, focus on answering the following research questions:

are different scalarizing functions perform different to each other under the same framework? 2. 2.

are surrogate models built on a given scalarizing function sensitive to the number of objectives ?

To answer the questions above, we embed different scalarizing functions into the framework of efficient global optimization (EGO) [31]. The algorithm uses the expected improvement criterion for updating the surrogates build on the scalarizing function. We then test the method with different scalarizing functions on several benchmark problems with different numbers of objectives. We analyze the accuracy and uncertainty provided by Gaussian process models on different scalarizing functions with the different number of objectives.

The rest of the article is organized as follows. In the next section, we provide a literature survey on using scalarizing functions with and without building surrogates on them. In Section 3, we give a brief introduction to different functions with their mathematical formulations and merits and demerits. In Section 4, numerical experiments are conducted to answer the research questions mentioned above. We finally conclude and mention the future research directions in Section 5.

2 Related work

Use of scalarizing functions has been a topic of interest in different communities. For instance, in Multiple Criteria Decision Making (MCDM) problems, several scalarizing functions have been used to incorporate a decision maker’s preferences [42, 21, 52, 43]. In evolutionary multiobjective optimization (EMO) algorithms especially decomposition based, utilization of scalarizing functions became popular in the last few years. Many studies and algorithms exist in the literature utilizing different functions for solving multi- and many-objectives optimization problems. Some of the well-known decomposition based EMO algorithms utilizing different scalarizing functions are nondominated-sorting genetic algorithm III (NSGA-III) [11, 27], multiobjective optimization based on decomposition (MOEA/D) [54] and its numerous versions [53] and reference vector guided evolutionary algorithm (RVEA) [4]. For more details on the working principle of decomposition based algorithms, see surveys on many-objective optimization evolutionary algorithms [24, 35].

In [26], two different scalarizing functions weighted sum and augmented Chebyshev was used adaptively in the solution process in the framework of MOEA/D by using a multi-grid scheme. The proposed idea was tested on a knapsack problem with four and six objectives and performed better than the original version of MOEA/D. The authors extended their work in [25], and modified the penalty boundary intersection (PBI) function in MOEA/D to handle different kinds of Pareto front. Two modified PBI functions called two-level PBI and quadratic PBI were tested in the MOEA/D framework and the algorithm with two new scalarizing functions performed better than the original version.

A detailed study on different scalarizing functions and their corresponding parameters was conducted in [13]. Instead of proposing a new algorithm, the authors showed the performance of the different scalarizing function in a simple (1+$\lambda$)-evolutionary algorithm [13] on bi-objective optimization problems. In [45], the performance of 15 different scalarizing functions was tested in MOEA/D [54] and MOMBI-II [18] algorithms. The authors used a tool called EVOCA [50] to tune the parameter values in different functions. Two algorithms MOEA/D and MOMBI-II [18] with different scalarizing functions were tested on Lame Supersphere test problems [14]. The authors found out that the performance of two algorithms depends on the choice of scalarizing functions.

In [16], a hyper-heuristic was used to rank different scalarizing functions with a measure called s-energy [17]. The proposed algorithm was tested on ZDT, DTLZ and WFG benchmark problems with 2-10 objectives and compared with MOMBI-II, MOEA/D and NSGA-III and found better results in most of the instances. In [28], two new scalarizing functions called the multiplicative scalarizing function (MSF) and penalty-based scalarizing function (PSF) was used in MOEA/D-DE [36]. The proposed scalarizing functions performed better than penalty boundary intersection, Chebyshev and weighted sum scalarizing functions.

The first algorithm in Bayesian multiobjective optimization method using a scalarizing function and building a surrogate on it was proposed in [33] and known as ParEGO (for Pareto based efficient global optimization). A Gaussian process model [49] was used as a surrogate of the scalarizing function. In the algorithm, a set of reference vectors was uniformly generated in the objective space using simplex lattice-design method [8]. In each iteration, the algorithm randomly selected a reference vector among a set of vectors, which was then used in the Chebyshev function to build the surrogate. The algorithm was compared with NSGA-II [12] and performed significantly better.

In a recent study in [48], three scalarizing functions, hypervolume improvement, dominance ranking, and minimum signed distance were proposed in solving expensive MOP. These scalarizing functions do not use reference vectors and preserve the dominance relationship. The authors also used the framework of EGO and compared with SMS-EGO [47] and ParEGO (i.e. with Chebyshev scalarizing function). The results outperformed the ParEGO algorithm and performed similarly to SMS-EGO in some cases. However, the reason for the better performance of their approach comparison to ParEGO was not mentioned.

3 Scalarizing functions

In this section, we summarize 15 different scalarizing functions studies in this work with their merits and demerits. All these functions have already been explained in details in the literature [7, 42]. Therefore, we provide a brief summary of these functions. Before going into description of scalarizing function, we define a multiobjective optimization problem (MOP) as:

[TABLE]

with $k(\geq 2)$ objective functions $f_{i}(x)$: S$\rightarrow\Re^{n}$. The vector of objective function values is denoted by $f(x)=(f_{1}(x),\ldots,f_{k}(x))^{T}$. The (nonempty) feasible space $S$ is a subset of the decision space $\Re^{n}$ and consists of decision vectors $x=(x_{1},\ldots,x_{n})^{T}$ that satisfy all the constraints. Different scalarizing functions are defined as follows:

Weighted sum (WS): The weighted sum combines different objectives linearly and has been widely used [42]. It converts the a MOP into a single-objective optimization problems as:

[TABLE]

One major limitation in using WS is that it cannot find solutions in non-convex parts of the Pareto front. 2. 2.

Exponential weighted criterion (EWC): EWC was first used in [2] to overcome the limitations of WS:

[TABLE]

As mentioned in [7], the performance of the function depends on the value of $p$ and usually a large value of $p$ is needed. 3. 3.

Weighted power (WPO): The weighted power defined below can find solutions when the Pareto front is not convex. However, it also depends on the parameter $p$ as in EWC.

[TABLE] 4. 4.

Weighted norm (WN): The weighted norm (or $L_{p}$ metric) is a generalized form of the weighted sum:

[TABLE]

The weighted norm has been used in MOEA/D [54]. Likewise EWC and WPO, its performance also depends on the choice of the $p$ parameter. 5. 5.

Weighted product (WPR): The weighted power is defined as:

[TABLE]

It is also called as product of powers [7] and can find solutions in non-convex parts of the Pareto front. However, like EWC, WPO and WN, its performance depends on the value of $p$. 6. 6.

Chebyshev function (TCH): The Chebyshev function can be derived from the WN function with $p=\infty$. This function has been used in many EMO algorithms such as MOEA/D [54] and its versions [53]:

[TABLE]

where $z_{i}^{*}$ is the ideal or utopian objective vector. In this work, we normalize the expensive objective function values in the range [0,1] before building the model. Therefore, $z_{i}^{*}$ is a vector of zeros. 7. 7.

Augmented Chebyshev (ATCH): In [52], it was suggested that weakly Pareto optimal solutions can be avoided by adding an augmented term to TCH:

[TABLE]

Moreover, this function was the first to be used as a surrogate in [33]. 8. 8.

Modified Chebyshev (MTCH): A slightly modified form of MTCH was used in [32]:

[TABLE]

As mentioned in [42] that main different between ATCH and MTCH is the slope to avoid weakly Pareto optimal solutions. 9. 9.

Penalty boundary intersection (PBI): The PBI function was first used in MOEA/D and used as the selection criterion for balancing convergence and diversity:

[TABLE]

where $d_{1}=|\mathbf{f}\cdot\frac{\mathbf{w}}{\|\mathbf{w}\|}|$ and $d_{2}=\|\mathbf{f}-d_{1}\frac{\mathbf{w}}{\|\mathbf{w}\|}\|$. The PBI function has been widely used in EMO algorithms [53]. However, as shown in [22] its performance is effected by the $\theta$ parameter. 10. 10.

Inverted penalty boundary intersection (IPBI): To enhance the diversity of solutions, IPBI was proposed in [51]:

[TABLE]

where $d_{1}=|\mathbf{f}^{{}^{\prime}}\cdot\frac{\mathbf{w}}{\|\mathbf{w}\|}|$ and $d_{2}=\|\mathbf{f}^{{}^{\prime}}-d_{1}\frac{\mathbf{w}}{\|\mathbf{w}\|}\|$ and $\mathbf{f}^{{}^{\prime}}=\mathbf{z}^{nadir}-\mathbf{f}$. In this work, we considered the vector of worst expensive objective function values as $\mathbf{z^{n}adir}$. 11. 11.

Quadratic PBI (QPBI): Recently, an enhanced version of PBI function was proposed in [25]:

[TABLE]

where $d_{1}$ and $d_{2}$ are same as in PBI and $d^{*}$ is an adaptive parameter and defined as:

[TABLE]

where $\alpha$ is a pre-defined parameter and $H$ is a parameter used in generating the reference (or weight vectors) in decomposition based EMO algorithms. For more details about these parameters, see [54, 4]. 12. 12.

Angle penalized distance (APD): The APD function is a recently proposed scalarizing function and used as the selection criterion in reference vector guided many-objective evolutionary algorithm (RVEA) [4]:

[TABLE]

where $P(\theta)=k(\frac{FE}{FE^{max}})^{\alpha}\frac{\theta}{\gamma}$, $FE$ is the number of expensive function evaluations at the current iteration and $FE^{max}$ is the maximum number of expensive function evaluations. The angle between an objective vector $\mathbf{f}$ and the reference vector to which it is assigned is represented by $\theta$ and the minimum of all angles between a reference vector selected and other reference vectors is represented by $\gamma$. This function adaptively balances the convergence and diversity based on the maximum number of function evaluations. The performance of the function depends on the value of $\alpha$ and $FE^{max}$. 13. 13.

Hypervolume improvement (HypI): The HypI is a recently proposed scalarizing function in a surrogate-based algorithm [48]. In contrast to its name, this scalarizing function is designed in such a way that it uses the Pareto dominance. Given a set of solutions $X$, a nondominated sorting is performed to find fronts of different ranks as in [12]. Let the different fronts of ranks $1,2,\ldots$ are denoted by $P_{1},P_{2},\ldots$ and the hypervolume of a front $P_{k}$ given a reference point $r$ is denoted by $H(P_{k})$. Then hypervolume improvement a solution $x\in X$ belongs to the front $P_{k}$ is then given by:

[TABLE]

In this way, the Pareto dominance is preserved when calculating the hypervolume contributions of solutions. For problems with many-objectives, where all solutions are non-dominated i.e. belong to only one front, this scalarizing function might not be suitable. In other words, if all solutions are nondominated, their HypI values will be similar and the algorithm embedding the model will not be able to solve problems with many-objectives. Moreover, the performance of the function can be sensitive to the reference point in calculating the hypervolume. 14. 14.

Dominance ranking (DomRank): This function is also recently proposed in [48] and assigns fitness values based on the ranks of different solutions as done in the MOGA algorithm [15]. In [15], a solution is assigned a rank as: rank(x)= 1+$p$, where $p$ is number of solutions dominating x. Similarly, in [48], given a set of solutions (with expensive evaluations) $X$ the fitness of a solution $x$ is:

[TABLE]

For instance, the rank of a solution $x^{\prime}$ dominated by all other solutions would be $rank(x^{\prime})=1+|X|-1$, and therefore, the fitness of solution would be $g(x^{\prime})=1-\frac{1+|X|-1-1}{|X|-1}=0$. Similarly, the rank of solutions belong the first front will be 1. This function is maximized to find samples for training the surrogates. Similar to HypI, this function also uses Pareto dominance and might not be suitable for problems with many-objectives. 15. 15.

Minimum Signed Distance (MSD) [48]: The MSD function is proposed in [48] and defined as:

[TABLE]

where $x^{\prime}$ are the solutions belong the first front (or rank one solutions) and $d(x^{\prime},x)=\sum_{i=1}^{k}f_{i}(x^{\prime})-f_{i}(x)$. For instance, if a solution a dominates another solution b, then $g(a,x^{\prime})>g(b,x^{\prime})$.

4 Numerical experiments

This section provides a comparison of different functions and analysis when building surrogates on them. As mentioned we used the framework of EGO [31] and built Gaussian process as surrogates on the scalarizing functions. The expected improvement criterion in EGO is maximized with a genetic algorithm for selecting samples when updating the surrogates.

4.1 Performance of different scalarizing functions

To compare the performance of different scalarizing functions, we used DTLZ [9] and WFG [20] problems with 2, 3, 5 and 10 number of objectives. In DTLZ suite, the number of variables is kept to $k$+5-1, where $k$ is the number of objectives. For the WFG suite, the number of variables is defined by position ($d$) and distance ($l$) parameters. We used the number of parameters as suggested in [20]. For two objectives, $d$ is set to four and $2\times(k-1)$ for the rest of the objectives. The distance parameter is set to four for all objectives. To be summarized, the number of variables ($n$) for DTLZ and WFG suites is given in Table 1:

There are several parameters in different scalarizing functions and we used the recommended values from the respective articles. Different parameter values used are provided in Table 2:

For WFG problems, the reference point $r$ in HyPI is used as $2\times z^{nadir}$. The vector $z^{nadir}$ contains the maximum objective function values in the Pareto front of the given problem. We ran 21 independent runs for each scalarizing function with 300 maximum number of expensive function evaluations. We show the performance of different functions with inverted generational distance (IGD) and hypervolume. The mean IGD values with standard deviation (in parentheses) for DTLZ problems are provided in Table 4 and other results including IGD and hypervolume on WFG are in the supplementary material111available from https://github.com/tichugh/SCF. To compare the results of different functions, we used the Wilcoxon rank sum test and used the Gaussian (or RBF) kernel when using the Gaussian process model.

In tables, the values statistically similar to the best value are in bold, the best value is in the circle and the worst value is underlined. As can be seen, the performance of different functions varies with different problems. Moreover, we provide in how many instances, a particular scalarizing function has the best IGD value in Figure 2. The total number of instances tested were (7 + 9) $\times$ 4 = 64 i.e. (number of DTLZ problems + number of WFG problems) $\times$ number of cases with respect to objectives.

As can be seen from tables and the figure, three functions PBI, HypI, and APD outperformed other functions. Moreover, WS, ATCH, and DomRank never had the best IGD value in any of 64 instances. Note that, all these three functions have statistically similar values to the best IGD value in many cases as can be seen from tables. Another interesting observation from the results is that three different scalarizing functions ATCH, DomRank and MSD used in the literature for building surrogates on them were not in the top list. On the other hand, scalarizing functions used as the selection criterion in EMO algorithms e.g. APD and PBI performed significantly better. These results indicate that it is not straightforward to select a particular scalarizing function. However, one can analyze the fitness landscape of different functions after building surrogates. Moreover, this landscape can be sensitive to the number of objectives and the details are provided in the next subsection.

4.2 Sensitivity towards number of objectives

When using scalarizing functions for building surrogates on them, the number of objectives can play a crucial role in their performance. Before using the surrogate built on the scalarizing function, we can analyze and visualize the fitness landscape with respect to the accuracy and uncertainty of approximations (or predictions). We start by showing the bar graph of the number of wins with respect to 2, 3, 5 and 10 objectives in Figure 3.

As can be seen in the figure that the performance of different functions varies with the number of objectives. The reason is that different formulations of scalarizing functions provide a different landscape. For a better illustration, an example of fitness landscapes of DTLZ2 with two and 10 objectives for different functions are shown in Figures 4 and 5, respectively. In figures, we used a training data set (shown as circles in figures) and built the surrogates on scalarizing functions. Then we generated a testing data set and used the model to get the predicted values (solid lines in figures) and uncertainties of the predicted values (shaded region in figures). The dotted line represents the actual scalarizing function values (denoted by g in figures) on the testing data set. For visualization purposes, we could show only one of the variables values.

As can be seen, for the same training data set (and therefore same objective function values), different scalarizing functions give different values. Therefore, the surrogate model to be built strongly depends on the g function values. For instance, in many-objective case, most of the solutions were nondominated and functions which use the property of Pareto dominance produce similar g values, as can be seen in the plots of HypI, DomRank, and MSD in 10 objective cases. As the g function values are the same for most of the solutions, the algorithm (EGO in this work) for finding a sample to update the surrogate by optimizing an infill criterion (EI in this work) will not be able to enhance the accuracy of the surrogates This is also the reason that HypI performed very well in a low number of objectives and its performance deteriorated with an increase in the number of objectives. On the other hand, functions like APD and PBI are not very sensitive to the number of objectives as they do not use the property of Pareto dominance and suitable for a large number of objectives when building surrogates on them. In addition to the sensitivity towards many objectives, some functions, EWC and WS did not perform well in most of the problems. In EWC, the g function landscape is very rugged [40] as can be seen in figures and WS has the problem of finding solutions in the non-convex parts of the Pareto front. These results show that one needs to see or analyze the landscape of scalarizing functions for the given training data set before using them in the optimization. We observed the similar behavior on other problems and landscapes on other problems with different objectives are provided in the supplementary material.

5 Conclusions and future research directions

This work focused on reviewing, analyzing and comparing different scalarizing functions in Bayesian multiobjective optimization for solving (computationally) expensive multi- and many-objective optimization problems. We provided an overview of different functions with their merits and demerits. We built the surrogates on the different functions in the framework of Bayesian multiobjective optimization. We compared different functions with the different number of objectives by using IGD and hypervolume as the performance metrics. The results clearly showed that some functions outperformed others in many cases and some did not work in most of the cases. We then analyzed the fitness landscape of different functions with respect to the number of objectives. We found out that some of the functions are very sensitive to the number of objectives and an analysis of the landscape of different functions might be helpful in selecting functions when building surrogates on them and using them in the optimization for solving expensive MOPs.

In this work, we did not compare the scalarizing functions with other non-Bayesian multiobjective optimization algorithms. Therefore, comparison with other algorithms are topics for future research directions. Moreover, working on an adaptive strategy in using and selecting scalarizing functions and their corresponding parameters are also in the list of future research.

Acknowledgement

This research is supported by the Natural Environment Research Council [NE/P017436/1]