1An Analysis of Quality Indicators Using
Approximated Optimal Distributions in a
Three-dimensional Objective Space
Ryoji Tanabe, Member, IEEE,and Hisao Ishibuchi, Fellow, IEEE

Abstract—Although quality indicators play a crucial role inbenchmarking evolutionary multi-objective optimization algorithms, their properties are still unclear. One promising approachfor understanding quality indicators is the use of the optimal
distribution of objective vectors that optimizes each quality indicator. However, it is difficult to obtain the optimal distribution foreach quality indicator, especially when its theoretical property is
unknown. Thus, optimal distributions for most quality indicators
have not been well investigated. To address these issues, first, we
propose a problem formulation of finding the optimal distribution
for each quality indicator on an arbitrary Pareto front. Then, we
approximate the optimal distributions for nine quality indicators
using the proposed problem formulation. We analyze the nine
quality indicators using their approximated optimal distributions
on eight types of Pareto fronts of three-objective problems.
Our analysis demonstrates that uniformly-distributed objective
vectors over the entire Pareto front are not optimal in manycases. Each quality indicator has its own optimal distribution foreach Pareto front. We also examine the consistency among thenine quality indicators.
Index Terms—Evolutionary multi-objective optimization, quality indicators, optimal distributions of objective vectors

I. I NTRODUCTIONM

ULTI-OBJECTIVE optimization problems (MOPs) appear in real-world applications. Since no solution cansimultaneously optimize multiple objective functions in general, the goal of MOPs is usually to find a Pareto optimal
solution preferred by a decision maker [1]. When no information about the decision maker’s preference is available, an “aposteriori” decision making is performed. The decision makerselects a final solution from a non-dominated solution set
that approximates the Pareto front in the objective space. An
evolutionary multi-objective optimization algorithm (EMOA)
is frequently used for the “a posteriori” decision making [2].
In multi-objective optimization, a vector that consists of all
objective values of a solution is referred to as an objective
vector. In this paper, we are interested in the distributionof objective vectors in the objective space, rather than the
distribution of solutions in the solution space. Two terms“objective vector” and “solution” are used synonymously insome previous studies (e.g., [3]–[5]). For the sake of clarity,
R. Tanabe and H. Ishibuchi are with Shenzhen Key Laboratory of Computational Intelligence, University Key Laboratory of Evolving Intelligent Systemsof Guangdong Province, Department of Computer Science and Engineering,Southern University of Science and Technology, Shenzhen 518055, China.
e-mail: (rt.ryoji.tanabe@gmail.com, hisao@sustech.edu.cn). (Corresponding
author: Hisao Ishibuchi)

we distinguish the two terms. That is, solutions in the objective
space are always referred to as objective vectors in thispaper. We also assume that the number of objective vectors isbounded by the population size µ.
Most quality indicators evaluate the quality of an objective
vector set in terms of at least one of the three criteria (theconvergence, the uniformity, and the spread) [5]–[8]. Althougha number of quality indicators have been proposed in the literature, we focus on unary quality indicators that map µ objective
vectors to a real number. Representative quality indicators
include hypervolume (HV) [9], the additive -indicator (I+ )[7], and inverted generational distance (IGD) [10]. Since theperformance of EMOAs is discussed based on the quality ofobtained objective vector sets, quality indicators play a central
role in benchmarking EMOAs. Although this paper does notaddress the cardinality of an objective vector set, it does notmean that the cardinality is less important. As pointed out in[5], the cardinality is one of important aspects when the size of
the Pareto front is relatively small (e.g., combinatorial MOPs[11]). Since we address only non-dominated objective vectors
on the Pareto front of a continuous MOPs in this paper, we
do not analyze the cardinality-based quality indicators.
One of the critical issues in quality indicators is that
their properties are not always obvious. Thus, some quality
indicators are likely to provide misleading information about
the quality of objective vector sets. HV is maximized by nonuniform objective vectors when the Pareto front is nonlinear[12]. Although the original generalized spread (∆) for biobjective problems [13] can evaluate the uniformity of given
objective vectors, its extended version for problems with morethan two objectives [14] overestimates the uniformity of nonuniform objective vectors in some cases [5], [15].
Since understanding the properties of quality indicators
is necessary for the development of EMOAs in the right
direction, some analysis studies have been presented in theliterature. One of the most promising approaches is the use of
the optimal distribution of objective vectors that optimizes agiven quality indicator. A distribution of µ objective vectors
is said to be the optimal µ-distribution if it optimizes aquality indicator [16]. The optimal µ-distribution can explain
which distribution of objective vectors is highly evaluated byeach quality indicator. For example, the optimal µ-distribution
for IGD does not contain objective vectors on the edges ofthe linear Pareto front in some cases [3]. Thus, IGD may
overestimate an objective vector set with a poor spread. The
optimal µ-distribution helps an indirect analysis of each quality

2

indicator in this manner.
Although the optimal µ-distributions for HV and IGD
can be obtained on the linear Pareto front of two objective
problems [3], [12], it is difficult to obtain the optimal µdistributions for other quality indicators in an exact manner,
especially problems with more than two objectives. Thus,
approximated optimal µ-distributions are used for the analysis
of quality indicators. In general, an approximated optimal µdistribution is obtained in an analytical or empirical manner[3], [12], [17]–[19]. However, existing analytical and empirical approaches have some limitations. Analytical approaches
utilize some specific properties of a quality indicator (e.g., thedecomposability of the HV calculation [18]). For this reason,
analytical approaches can be applied only to quality indicators
whose theoretical properties are clear (e.g., HV for two- andthree-objective problems [12], [18] and I+ only for twoobjective problems [19]). Empirical approaches translate theproblem of finding the optimal µ-distribution on the Pareto
front into a black-box single-objective continuous optimization
problem. Then, the optimal µ-distribution is approximatedusing a derivative-free optimizer. For example, the optimalµ-distribution for R2 was approximated by CMA-ES [20] in[17]. Although empirical approaches are more flexible thananalytical approaches, empirical approaches have been applied
to only a few quality indicators, including HV [12] and R2[17]. In addition, most previous studies addressed only twoobjective problems. This is because it is not obvious how toformulate the problem of finding the optimal µ-distribution
on a multi-objective problem with more than two objectives.
For these reasons, the (approximated) optimal µ-distributions
for most quality indicators on problems with more than two
objectives have not been well analyzed in the literature.
To address these issues, first, we propose a problem formulation of finding the optimal µ-distribution for each quality indicator on a Pareto front in an arbitrary-dimensional objective
space. Then, we approximate the optimal µ-distribution foreach of the following nine quality indicators on a Pareto front
in a three-dimensional objective space: HV, IGD, IGD plus
(IGD+ ) [21], R2, new R2 (NR2) [22], I+ , s-energy (SE) [23],∆, and pure diversity (PD) [24]. We analyze the nine quality
indicators using their approximated optimal µ-distributions oneight types of Pareto fronts.
Our main contributions in this paper are as follows:
1) We propose a general problem formulation of finding the
optimal µ-distribution. In contrast to the existing empirical approaches, the proposed formulation can handle
any number of objectives. The proposed formulation isapplicable to any types of the Pareto front (under thecondition that its front shape function is given). Note thatthe formulation proposed in [12] and our formulation
can be categorized into set-based optimization [25].
Details are discussed in Section IV.
2) We approximate the optimal µ-distributions for the nine
quality indicators on the Pareto fronts of eight problems
with three objectives. This is the first study to present
the approximated optimal µ-distributions for NR2, ∆,SE, PD, R2, and I+ for three-objective problems.

3) We analyze the nine quality indicators from the viewpoint of the approximated optimal µ-distributions. We
provide insightful information about the nine quality indicators. For example, PD may overestimate an objective
vector set with a small dissimilarity. Our observations
are summarized in Section VII.
The rest of this paper is organized as follows. Section
II provides some preliminaries. Section III describes related
studies. Section IV proposes a problem formulation to search
for the optimal distribution of objective vectors. Section V describes the setting of our computational experiments. Section
VI shows analysis results. Section VII concludes this paper.
II. P RELIMINARIES
A. Definition of multi-objective optimization problems (MOPs)
We suppose a multi-objective minimization problem in thispaper. A continuous MOP is to find a solution x ∈ S thatminimizes a given objective function vector f : Rn → Rm .
Here, S ⊆ Rn is the n-dimensional solution space, and Rm is
the m-dimensional objective space. Thus, n is the number ofdecision variables, and m is the number of objective functions.
A solution x1 is said to dominate x2 iff fi (x1 ) ≤ fi (x2 )for all i ∈ {1, ..., m} and fi (x1 ) < fi (x2 ) for at least one
index i. In addition, x1 is said to weakly dominate x2 iff
fi (x1 ) ≤ fi (x2 ) for all i ∈ {1, ..., m}. If x∗ is not dominatedby any other solutions in S, x∗ is a Pareto optimal solution.
The set of all x∗ is the Pareto optimal solution set. The setof all f (x∗ ) is the Pareto front. The goal of MOPs for the “aposteriori” decision making is to find a non-dominated solution
set that approximates the Pareto front in the objective space.
B. Quality indicators
For the sake of simplicity,
vectorT we denote the objectivef (x) = f1 (x), ..., fm (x)as a = (a1 , ..., am )T , whereai = fi (x) for i ∈ {1, ..., m}. Let A = {a1 , ..., aµ } be a setof µ objective vectors, where µ is the population size.
Let Ω be a set of all possible objective vector sets. A unary
quality indicator I : Ω → R evaluates the quality of A interms of at least one of the convergence, the uniformity, andthe spread. I is said to be Pareto-compliant iff the rankingof all objective vector sets in Ω by I is consistent with thePareto dominance relation (for details, see [26]). In this paper,
we define the convergence, the uniformity, and the spread as
follows. The convergence of A means how close objective
vectors in A are to the Pareto front. The uniformity of A
means how uniform the distribution of objective vectors in A
is. The spread of A means how well objective vectors in A
cover the entire Pareto front. As pointed out in [5], the spread
is not equal to the extensity, which represents the length of
the boundaries of A. A combination of the uniformity andthe spread is referred to as “diversity” in the evolutionary
multi-objective optimization (EMO) community. Fig. S.1 inthe supplementary file shows examples of distributions ofobjective vectors in the three cases on the linear Pareto front
with m = 2.
Below, we explain the following nine quality indicators
used in this paper: HV, IGD, IGD+ , R2, NR2, I+ , SE,3∆, and PD. We do not consider convergence-based quality
indicators, such as generational distance (GD) [27], which
evaluate the quality of an objective vector set A in terms onlyof the convergence. This is because the optimal µ-distribution
for such a quality indicator is obvious. A is optimal for
convergence-based quality indicators if all elements in A areon the Pareto front regardless of their distribution. For thesame reason, we do not consider cardinality-based quality
indicators, such as overall non-dominated vector generation(ONVG) [28], which evaluate the quality of A based on thenumber of non-dominated solutions. For example, all of A1 ,A2 , and A3 in Fig. S.1 are optimal for both GD and ONVG.
Below, R is a set of m-dimensional reference vectorsthat approximate the Pareto front in the objective space. The
reference point q = (q1 , ..., qm )T ∈ Rm is used for the HV
and NR2 calculations. It should be noted that q 6∈ R.
1) Hypervolume (HV): Since HV [9] is the only Paretocompliant quality indicator known so far [29], HV is one of
the most popular quality indicators. The HV value of A is thevolume of the area that is dominated by objective vectors in
A and bounded by the reference point q as follows:!
[
HV(A) = volume[a1 , q1 ] × ... × [am , qm ] , (1)a∈Awhere the function “volume” in (1) is the Lebesgue measure.
A large HV value indicates that A approximates the Pareto
front well in terms of both convergence and diversity.
2) Inverted generational distance (IGD): IGD [10] evaluates the quality of A in terms of both convergence and diversity. IGD measures the average distance from each reference
vector r in R to its nearest objective vector a in A as follows:


n
o
1 X
(2)
min d(a, r)  ,IGD(A) =
|R|
a∈A
r ∈R
where d(a, r) in (2) is the distance between a and r.
IGD
pPmuses the following Euclidean distance: dIGD (a, r) =2
i=1 (ai − ri ) .
3) IGD plus (IGD+ ): Since IGD is not Pareto compliant,
IGD may incorrectly evaluate the quality of A that contains
non-converged objective vectors [21], [30]. IGD+ [21] addresses the issue in IGD. While IGD is Pareto non-compliant,
IGD+ is weakly Pareto compliant. The IGD+ value of A
is the average distance from each reference vector r to its
nearest objective vector in the dominated region by A [31].
Although the IGD+ value of A is calculated using (2), d is
replacedby the following distance function: dIGD+ (a, r) =q
2
Pm
i=1 max{ai − ri , 0} . It should be noted that small IGD
and IGD+ values are preferable whereas large HV values are
preferable.
In addition to IGD+ , the averaged Hausdorff distance indicator (∆p ) [30] addresses the issue in IGD. When all elements
in A are on the Pareto front and R has a large number of
uniformly distributed reference vectors over the entire Pareto
front (as in our computational experiments in this paper), the∆p value of A is equal to the IGD value of A in general.
Actually, the IGD and ∆p values are the same for all of theobtained solution sets in this paper. For this reason, we do
not report experimental results for ∆p . Of course, the abovementioned relationship between ∆p and IGD does not always
hold depending on A and R.
4) The additive -indicator (I+ ): Here, we describe theunary version of I+ [7]. The I+ value of two objective
vectors a and r is given as follows:vecI+(a, r) =maxi∈{1,...,m}{ai − ri },

(3)

vec
where the I+(a, r) value is the minimal shift such that avecweakly dominates r. It should be noted that I+(a, r) 6=
vec
I+ (r, a) in (3). I+ evaluates how well A approximates thereference vector set R as follows:

n
o
vec
I+ (A) = max min I+(a, r) ,(4)r ∈R a∈A
where the I+ value of A in (4) is the minimal shift such thateach reference objective vector in R is weakly dominated byat least one objective vector in A. A small I+ value indicatesthat the corresponding A is a good approximation of the Pareto
front in terms of both convergence and diversity.
5) R2: Although three R indicators are proposed in [32],the following R2 indicator is the most widely used one:noX1min g tch (a, w) ,(5)
R2(A) =
|W |
a∈A
w ∈W

g tch (a, w) = max
wi |ai − zi∗ | ,(6)i∈{1,...,m}where W in (5) is a set of weight vectors. For each w =(w1 , ..., wm )T ∈ W , 0 ≤ wi ≤ 1 for i ∈ {1, ..., m} andPmtch: Rm → R is the weighted
i=1 wi = 1. In (6), g∗∗ T) is the ideal
Tchebycheff function, and z = (z1∗ , ..., zm∗∗
point. The i-th element zi of z is the minimum value ofthe i-th objective function over the Pareto front. According to[5], R2 with a tuned z ∗ works well on various problems with
m = 2. A small R2 value indicates that A approximates thePareto front well in terms of both convergence and diversity.
6) New R2 indicator (NR2): Although the Pareto compliance of HV is attractive, the computation time for the HV
calculation increases exponentially with the number of objectives m. To address this issue, the method of approximating
the HV value using R2 is proposed in [33]. The following
NR2 [22] is an improved version of R2 to obtain a better
approximation of the HV value:nomX1
NR2(A) =min g mtch (a, w) ,
(7)
|W |
a∈A
w ∈W


|qi − ai |g mtch (a, w) = max,(8)wii∈{1,...,m}where g mtch in (8) is a modified version of g tch . If wi = 0
in (8), wi was set to 10−6 to avoid division by zero. Similar
to HV, NR2 uses the reference point q. A large NR2 value
indicates that the corresponding A has a good HV value. Since
NR2 was designed to approximate the HV value, it is expectedthat their optimal µ-distributions are similar to each other.
However, such an expectation has not been verified in theliterature. To address this issue, we examine the approximatedoptimal µ-distribution for NR2 in this paper.

47) s-energy (SE): Although SE [23] was not originallyproposed in the field of EMO, SE is used as a quality indicatorin some studies (e.g., [34]). SE is given as follows:
X
X
SE(A) =
ka − bk−s(9)2 ,a∈A b∈A\{a}where s is usually set to m−1, and kakp is the Lp norm of a.
If a = b in (9), we set the SE value of A to an infinitely large
value to avoid division by zero. A small SE value indicatesthat A has a good uniformity and a good spread.
8) Generalized Spread (∆): The original ∆ indicator [13]can handle only the two-dimensional objective space. In thispaper, we use an extension of ∆ with respect to m presented in[14]. ∆ evaluates the quality of A in terms of both uniformity
and spread as follows:P
dext + a∈A |dist(a, A) − davg |,(10)∆(A) =
dext + davg (|A| − m)dist(u, V ) =min ku − vk2 ,(11)
v ∈V ,v 6=umX
dext =
dist(r ext(12)
i , A),
i=1

d

avg

=

1 X
dist(a, A),
|A|
a∈A

(13)

where the function dist(u, V ) in (11) returns the Euclidean
distance from u to its nearest vector in the set V . In (12),r ext
is the i-th extreme objective vector with the maximum
ivalue of the i-th objective function in R (i ∈ {1, ..., m}). In(10), dext aims to evaluate the spread of A, while the otherparts aim to evaluate the uniformity of A. A small ∆ value
indicates that A has good spread and uniformity.
9) Pure diversity (PD): PD [24] is an extended version of
the Solow-Polasky diversity indicator [35] for multi-objective
optimization. The PD value of A is based on the dissimilarityof each objective vector a to A as follows:noPD(A) = max PD A \ {a} + d a, A \ {a} , (14)a∈A d(u, V ) = min ku − vkp ,
(15)
v ∈V
where the recommended value of p in the Lp norm in (15) is0.1. Most diversity-based quality indicators (including SE and∆) evaluate the uniformity of A. In contrast, PD evaluates thedissimilarity between objective vectors in A.

The other is optimal µ-distribution based approaches. As
explained in Section I, the optimal µ-distribution provides
insightful information about the properties of each quality
indicator. One main drawback of the optimal µ-distribution
based approaches is that it is difficult to find the optimal µdistributions for quality indicators.
A. Ranking information based approaches
Okabe et al. [8] analyzed 15 quality indicators using objective vector sets artificially generated on the linear Pareto front
and objective vector sets found by EMOAs on the convex
Pareto front. The results show that some quality indicators
can generate misleading results. Jiang et al. [15] examined six
quality indicators using a pair of objective vector sets on thelinear, convex, and concave Pareto fronts. The results showthat I+ , IGD, and HV are consistent with each other on theconvex Pareto front, but IGD is inconsistent with HV on the
concave Pareto front. Ravber et al. [36] investigated 11 quality
indicators using a chess rating system. They used objective
vector sets found by EMOAs on some test problems. The
results show the robustness of HV, I+ , and IGD+ . However,
the results also show that the three quality indicators do notequally evaluate the quality of objective vector sets in termsof the convergence, the uniformity, and the spread.
Wessing and Naujoks [37] analyzed the correlation between
HV, R2, and I+ . Their analysis was based on objective vector
sets randomly generated on two- and three-objective problems.
The results based on Pearson’s correlation coefficient showthat HV is highly consistent with R2. Liefooghe and Derbel[38] examined the correlation between six quality indicators.
They used objective vector sets found by the random searchand NSGA-II [13] on test problems with m ∈ {2, 3}. The
results based on the Kendall rank correlation show that HV
is consistent with R2 and R3. In contrast, I+ and its multiplicative version are inconsistent with other quality indicators.
Ishibuchi et al. [39] compared IGD+ with IGD on various
objective vector sets found by NSGA-II on problems with up
to 10 objectives. The results show that IGD+ can evaluate
objective vector sets more accurately than IGD in terms of thePareto compliance. In addition to IGD+ and IGD, three quality
indicators (GD, I+ , and D1 [32]) were examined. Bezerra etal. [40] analyzed the relation between IGD+ and IGD on thelinear Pareto front. Various objective vector sets were usedin their analysis. The results show that IGD and IGD+ are
consistent with each other in most cases. The results also showthat IGD may favor objective vectors with a poor spread.

III. R ELATED WORKThis section describes related work. Analysis of quality
indicators can be categorized into two approaches. One isranking information based approaches. First, multiple objective vector sets are ranked by each quality indicator. Then,the properties of each quality indicator are discussed based on
which objective vector set is highly evaluated by the quality
indicator. The consistency of multiple quality indicators canalso be investigated in this manner. One main drawback of
the ranking information based approaches is that only relative
information can be obtained.

B. Optimal µ-distribution based approaches
It is difficult to obtain the optimal µ-distributions in anexact manner. Thus, the optimal µ-distributions are generally
approximated by analytical or empirical methods. For thedifferences between the two methods, see Section I.
Auger et al. [12] formulated finding the optimal µdistribution for HV on two-objective problems as a singleobjective continuous optimization problem. Below, we describe its problem formulation in a general manner. Our
description significantly differs from the original one.

5

0.5

0
0

0.5

(a) FDTLZ1

1

1

0.5

0.5

0
0

0.5

(b) FDTLZ2

1

0
0

0.5

1

(c) FZDT1Fig. 1: Distributions of five translated objective vectors. The x andy axes represent f1 (θ) and f2 (F (θ)), respectively.

Let A = {a1 , ..., aµ } be a set of µ objective vectors on the
Pareto front. Also, let θi be a real value that determines thefirst element (i.e., the value of the first objective f1 ) of the i-th
objective vector ai (i ∈ {1, ..., µ}) in A. All objective vectors
in the optimal µ-distributions must always be on the Pareto
front. For this restriction, the position of ai for f1 specifies
that for the second objective f2 . In [12], a front shape function
F : R → R is used to determine
the position of ai for f2 asTfollows: ai = θi , F (θi ) . For all i ∈ {1, ..., µ}, θi is in
the range [θmin , θmax ], where θmin and θmax are the lower
and upper bounds for θi . Here, θmin and θmax correspond tothe minimum and maximum values of f1 in the Pareto front,
respectively. In this manner, the distribution of µ objective
vectors in A is determined by a vector θ = (θ1 , ..., θµ )T .
Fig. 1 shows examples of five objective vectors. In Fig. 1,
the Pareto front is shown by
FDTLZ1 (θ) = 0.5 − θp
FDTLZ2 (θ) = 1 − θ2√
FZDT1 (θ) = 1 − θ(θ ∈ [0, 0.5]),(16)(θ ∈ [0, 1]),(17)

(θ ∈ [0, 1]),(18)

where FDTLZ1 , FDTLZ2 , and FZDT1 represent the Pareto
front shapes of DTLZ1 [41], DTLZ2 [41], and ZDT1 [42],respectively. In Fig. 1, µ = 5, θ = (0, 0.125, 0.25, 0.375, 0.5)Tfor FDTLZ1 , and θ = (0, 0.25, 0.5, 0.75, 1)T for FDTLZ2 andFZDT1 . As shown in Fig. 1, distributions of objective vectors
on various Pareto fronts can be examined by changing F .
In summary, finding A that minimizes a quality indicator I
is the same as finding θ that indirectly minimizes I as follows:noTminimize I(A) = I
θi , F (θi ) | i ∈ {1, ..., µ} , (19)where θ is a solution of the problem defined in (19). I thatneeds to be maximized (e.g., HV, NR2, and PD) can betranslated as −I without loss of generality. Any derivativefree optimizer can be used to minimize I in (19).
Auger et al. [12] proposed analytical and empirical methodsof approximating optimal µ-distributions for HV. Details of
the analytical method for each F can be found in the supplementary website of [12] (https://sop.tik.ee.ethz.ch/download/
supplementary/testproblems/). The proposed problem-specificlocal search algorithm in [12] for (19) uses variable-wise
perturbation. For each iteration, each θi is perturbed by a value
selected randomly from a normal distribution. Brockhoff etal. [17] applied CMA-ES [20] to (19) to find the optimal µdistributions for R2. Their results show that objective vectors

in the optimal µ-distributions for R2 are densely distributedin the center of the Pareto front. They also examined the contribution of each objective vector in the optimal µ-distribution
for R2 to HV (and that for HV to R2).
In contrast to the above-mentioned studies [12], [17], thefollowing studies do not use (19). Bringmann et al. [19]proposed the analytical method for finding the optimal µdistributions for I+ on the Pareto front with m = 2.
The proposed method explicitly exploits the properties ofI+ presented in [43]. The results show that the optimal µdistributions for I+ on most Pareto fronts with m = 2 are
uniform. Glasmachers [18] proposed a gradient-based analytical method for approximating the optimal µ-distributions for
HV on m ∈ {2, 3}. For m = 3, the proposed method in[18] decomposes the three-dimensional objective space into
cuboids. Then, the proposed method in [18] uses the gradientof HV based on the cuboids.
Ishibuchi et al. [4] used SMS-EMOA [44] to approximate
the optimal µ-distributions for HV on Pareto fronts of threeobjective problems. SMS-EMOA uses the HV contribution
in the steady-state selection. For each iteration, the worst
individual in terms of the HV contribution is removed fromthe last non-domination level. In [4], SMS-EMOA was applied
to the DTLZ and inverted-DTLZ problems whose numberof distance variables is zero. Then, the influence of thereference vector on the optimal µ-distributions for HV is
examined. Ishibuchi et al. [3] also investigated the optimal µdistributions for IGD approximated by a similar approach. The
HV contribution-based selection in SMS-EMOA was replacedby the IGD contribution-based selection. The optimal µdistributions for HV, IGD, and IGD+ approximated by SMSEMOA in a similar manner were analyzed in [31]. The results
show that the optimal µ-distributions for IGD are almost
uniform in all problems. The results also show that the optimalµ-distributions for HV and IGD+ are similar in all problems
when the reference vector for HV is set appropriately.
IV. P ROPOSED FORMULATION
Here, we introduce the proposed problem formulation offinding the optimal µ-distribution for quality indicators on
Pareto fronts for arbitrary number of objectives (m ≥ 2). As
reviewed in Subsection III-B, some methods of approximating
the optimal µ-distributions have been proposed in the literature. The analytical methods [12], [18], [19] can efficiently
approximate the optimal µ-distributions in a computationallycheap manner. Unfortunately, they can be applied to only
quality indicators whose theoretical properties are clear (i.e.,
HV and I+ ). It is questionable whether the empirical methods
using SMS-EMOA [3], [4], [31] can be applied to quality
indicators except for HV, IGD, and IGD+ . The empiricalmethods based on the problem formulation in (19) are promising. Although the previous studies [12], [17] addressed only
HV and R2, the approaches using (19) can be applied toany quality indicators in principle. However, the problem
formulation in (19) can handle only Pareto fronts for m = 2. In
contrast, the proposed problem formulation can handle Pareto
fronts for m ≥ 2. We do not claim that the degenerated version

6of our formulation for m = 2 is better than the formulation in(19). Note that our formulation for m = 2 is almost the sameas the formulation in (19).
The formulation proposed in [12] for m = 2 and ourformulation for m ≥ 2 are essentially one of the setbased optimization methods [25]. Test problems constructedby bottom-up approaches (e.g., DTLZ) have two types ofdecision variables [45]. One is position-related variables thatspecify the position of an objective vector on the Pareto
front. The other is distance-related variables that determine
the distance between the objective vector and the Pareto front.
Thus, any solutions with only position-related variables arealways on the Pareto front in the objective space. The optimalµ-distribution for a quality indicator can be approximated by a
set-based method that finds a set of µ solutions (not objective
vectors) with only position-related variables for minimizationof the quality indicator. While we here denote this approach as
the solution set-based approach, we denote the formulations
proposed in [12] and this paper as the objective vector setbased approach. The main advantage of the objective vector
set-based approach is its simplicity. While the solution setbased approach needs to handle both solution and objective
spaces, the objective vector set-based approach just needs tohandle only the objective space. We can focus only on theobjective space in the objective vector set-based approach.
This simplicity is beneficial for an analysis of the optimalµ-distributions as demonstrated in [12].
A. Proposed problem formulation
Below, we explain the proposed formulation for generalcases. Only a disconnected front shape function Fdisconnected
explained in Subsection IV-C requires an alternative formulation. For details, see Subsection S.1-D in the supplementary
file. Finding the optimal µ-distribution for a quality indicator I
on the Pareto front with m objectives is formulated as follows:
minimize I(A),(20)

subject to G(A) ≤ 0,
where G is a constraint function for constrained front shapefunctions (e.g., Fc-concave explained in Subsection IV-C). In
contrast to (19), G is added to (20). If a given front shapefunction has no constraint, G is eliminated from (20). An
objective vector set A = {a1 , ..., aµ } is a translated version ofa solution θ of this problem in a two-phase manner as follows:
A = F (B),(21)

B = translate(θ),

(22)

T

where θ = (θ1 , ..., θd ) is a d-dimensional vector, and d =µ(m − 1). Any single-objective black-box optimizer can beused to find θ that minimizes I. While the size of θ is µ in(19), that is d in (21). For all i ∈ {1, ..., d}, θi is in [0, 1].
The function “translate” in (22) translates a d-dimensional
vector θ into a set B = {b1 , ..., bµ } that consists ofµ m-dimensional vectors.
Pm For each i ∈ {1, ..., µ}, bi =(bi,1 , ..., bi,m )T , and
j=1 bi,j = 1. Although bi is an mdimensional vector, it is not an objective vector. Subsection
IV-B explains the translation operator in detail.

After translating θ into B in (22), B is further mapped tothe objective vector set A by a front shape function F in (21).
For each i ∈ {1, ..., µ}, an element bi in B corresponds to an
objective vector ai in A. All objective values are normalized
into [0, 1] in our study by using the ideal point z ∗ and thenadir Tnadir point z nadir = (z1nadir , ..., zm) of each F as follows:
∗
nadir
∗
ai := (ai − zi )/(zi− zi ) for each i ∈ {1, ..., m}.
Here, zinadir is the maximum value of the i-th objectivefunction of the Pareto front. We use the normalization inorder to examine the approximated optimal µ-distributions forthe quality indicators without worrying about the influence ofdifferently-scaled objective values. Subsection IV-C specifies
how to map B to A by the front shape function F .
B. Translation operator
We explain how the translation operator in (22) maps θ
to B. First, all d elements in θ are divided into µ vector
groups as (θ1 , ..., θ1+m−2 )T , ..., (θd−m+2 , ..., θd )T . For example, when µ = 2, m = 3, d = µ(m − 1) = 4, andθ = (0.1, 0.2, 0.3, 0.4)T , 4 elements in θ are divided asfollows: (0.1, 0.2)T and (0.3, 0.4)T . For the sake of simplicity,
we denote (θi , ..., θi+m−2 )T as y = (y1 , ..., ym−1 )T for eachi ∈ {1, 1(m − 1), 2(m − 1), ..., (d − m + 2)}.
Below, we describe how the (m − 1)-dimensional vector yis translated into an m-dimensional vector b in B. To obtain
all µ vectors in B, the same operation is repeatedly applied
to y for each i ∈ {1, 1(m − 1), 2(m − 1), ..., (d − m + 2)}.
Each element in b is obtained using the Jaszkiewicz’s methodof generating random weight vectors [46] as follows:√(23)
b1 = 1 − m−1 y1 ,
!
j−1
X
√ 
bj = 1 −bk 1 − m−1−j yj ,bm = 1 −

k=1
m−1
X

bk ,

k=1

where j ∈ {2, ..., m − 1} in (23). It should
Pmbe noted that b
obtained by (23) satisfies the condition
k=1 bk = 1. Although the Jaszkiewicz’s method was originally proposed forthe weight vector generation in decomposition-based EMOAs,
it can be used to translate y into b in a straightforward manner.
C. Front shape function F
The proposed formulation can be applied to any Pareto front
when its front shape function F is available. In other words,
the proposed formulation cannot be applied to the Pareto front
whose analytical expression (i.e., its front shape function)is unavailable such as WFG3 [47]. Although various front
shape functions were presented in the previous study [12],they are only for m = 2. We design front shape functions form ≥ 3 based on a method of generating uniformly distributed
reference vectors in representative test problems (e.g., DTLZ
and WFG) presented in [48].
We use eight front shape functions shown in Table I. The
name of the original problems and the shape of the Pareto
fronts are also shown in Table I. Flinear , Fconcave , and Fconvex

7

TABLE I: Properties of the eight front shape functions.
F

Original problem

Shape

Flinear
Fconcave
Fconvex

DTLZ1 [41]DTLZ2 [41]Convex DTLZ2 [49]LinearConcave
ConvexFi-linearFi-concave
Fi-convex

Inverted DTLZ1 [50]Inverted convex DTLZ2 [51]Inverted DTLZ2 [51]Inverted, linear
Inverted, concave
Inverted, convex

Fdisconnected
Fc-concave

DTLZ7 [41]C2-DTLZ2 [50]Disconnected, mixed
Disconnected, concave

are the most basic front shape functions. Fi-linear , Fi-concave ,and Fi-convex are inverted versions of Flinear , Fconvex , andFconcave , respectively. It should be noted that an inverted
version of the convex Pareto front is concave, and vice versa.
Fdisconnected has the disconnected and mixed Pareto front.
Some parts of the Pareto front of Fc-concave are infeasible dueto the constraint. We select Fc-concave to demonstrate that the
proposed formulation can handle the constrained Pareto front.
In our preliminary experiments, we used the degenerate front
shape function Fdegenerate whose original problem is DTLZ5.
Fdegenerate is not a surface but a curve since it is degenerate.
Because the shape of Fdegenerate is concave, the optimal µdistribution on Fdegenerate is similar to that for a two-objective
problem with Fconcave . For this reason, we omit Fdegenerate .
Due to the paper length limitation, we explain only Fconcavein this paper. We describe other front shape functions in
Section S.1 in the supplementary file. In Fconcave , the i-th
element bi of b in B is translated into the i-th element ai of
a in A on the concave Pareto front as follows:bi(24)
ai = ,tqPm2
where t =j=1 bj . This translation is based on the methodof generating reference vectors in DTLZ2 described in [48].
All objective vectors obtained by (24) P
satisfy the followingmrelation of the Pareto front of DTLZ2: i=1 (fi (x∗ ))2 = 1.
In addition to the eight front shape functions in Table I, we
believe that it is possible to design other front shape functions
using [48] as reference.
V. E XPERIMENTAL SETTINGS
This section describes experimental settings. We examine
the properties of the nine quality indicators using the eight
front shape functions in Table I. Although the number ofobjectives m can be arbitrarily specified in the proposed
formulation, m is set to three in order to visually discuss
the distributions of objective vectors. We investigate the optimal µ-distributions with µ = 10, 15, 21, 28, 36, 45 (d =20, 30, 42, 56, 72, 90, respectively). We select these µ values
so that we can visually discuss the distribution of objective
vectors. The distribution of objective vectors with a too largeµ value is unclear. We can also examine the influence of µ on
the optimal µ-distributions by using the various µ values.
In the proposed formulation in (20), a black-box optimizer
is necessary to find θ that minimizes a given quality indicator.

We use L-SHADE [52], which is an improved version ofdifferential evolution (DE) [53]. L-SHADE was the winnerof the IEEE CEC2014 competition on single-objective realparameter optimization. We used the Java source code of
L-SHADE provided by the authors of [52]. The default
parameter setting was used for L-SHADE. The maximum
number of function evaluations was set to 104 × d, and 31
independent runs were performed. For the constrained front
shape function Fc-concave , we use the death-penalty constraint
handling method. This is because objective vector sets are
just classified into feasible and infeasible groups as describedin Subsection S.1-E in the supplementary file. The fitness
value of an infeasible individual is an infinitely large value
so that all infeasible individuals cannot survive to the next
iteration. We compared L-SHADE to the analytical approachthat approximates the optimal µ-distributions for HV form = 2 [12] (see Subsection III-B). We confirmed that LSHADE can find the approximated optimal µ-distribution for
HV with acceptable quality for m = 2. For details, see Section
S.3 in the supplementary file.
We implemented the nine quality indicators ourselves, except for HV. The WFG algorithm [54] was used for the
HV calculation. We used the source code of the WFG algorithm provided by the authors of [54]. For PD, we carefullytranslated the Matlab source code provided by the authors of[24] into our implementation in order to speed up the PD
calculation. We confirmed that our version and the original
version output exactly the same PD value on all the front
shape functions used in our study. The source codes used inour experiments can be downloaded from the supplementary
website (https://sites.google.com/view/optmudist).
We set the reference vector q for HV and NR2 to(1.2, ..., 1.2)T according to [4]. We set z ∗ for R2 to (0, ..., 0)T .
As mentioned in Subsection IV-A, all objective values are in
the range [0, 1] for all front shape functions. For R2 and NR2,
we generated the weight vector set W using simplex-latticedesign [55]. For IGD, IGD+ , I+ , and ∆, we generated thereference vector set R by applying each front shape function F
to W as follows R = F (W ). That is, in order to generate R,
W is used in (21) instead of B (the resulting A is used as R).
For details, see Section S.2 in the supplementary file. The sizeof W and R was set to 1 035. For Fdisconnected and Fc-concave ,we set the size of R to 1 089 and 1 087, respectively. This isdue to the properties of Fdisconnected and Fc-concave .

VI. E XPERIMENTAL RESULTS
This section analyzes the nine quality indicators using theiroptimal µ-distributions approximated by the proposed formulation. Subsection VI-A discusses the optimal µ-distributions
for each quality indicator. Subsection VI-B examines the nine
quality indicators using the ranking information. Subsection
VI-C analyzes the generality of results for m = 3 with respect
to m. Subsection VI-D investigates a unary version of an M nary quality indicator.

8(a) AHV(b) AIGD(a) AHV(b) AIGD(a) AHV(b) AIGD(c) AIGD+(d) AR2(c) AIGD+(d) AR2(c) AIGD+(d) AR2(e) ANR2(f) AI+(e) ANR2(f) AI+(e) ANR2(f) AI+(g) ASE(h) A∆(g) ASE(h) A∆(g) ASE(h) A∆(i) APD(j) ASLD(i) APD(j) ASLD(i) APD(j) ASLD

Fig. 2: Results on Flinear .

Fig. 3: Results on Fconcave .

A. Approximated optimal µ-distributions for the nine quality
indicators
Figs. 2–9 show approximated optimal µ-distributions with
µ = 21 obtained by the proposed approach for the nine quality
indicators on the eight Pareto fronts. In Figs. 2–9, x, y, andz axes represent f1 , f2 , and f3 , respectively. Figs. 2–9 show
the best objective vector set found by L-SHADE among 31
runs for each quality indicator on each front shape function F .
Below, we denote the best µ-distribution for a quality indicator
I as AI (e.g., AHV ). Figs. S.3–S.50 in the supplementary fileshow approximated optimal µ-distributions with the other µvalues. All reference vectors in R are shown in each figure.
For the sake of comparison, Figs. 2–7 (j) show a setof µ objective vectors generated by simplex-lattice design,
denoted as ASLD . We generated ASLD using the same method
as for generating the reference vector set R described in
Section V. Since the number of objective vectors generated bysimplex-lattice design on Fdisconnected and Fc-concave cannot
be specified arbitrarily due to their characteristics, we do notshow ASLD in Figs. 8 and 9. ASLD can be viewed as an“ideal” objective vector set found by decomposition-based
EMOAs (e.g., NSGA-III [49]). Although the distribution of
objective vectors in ASLD is uniform, ASLD is not optimal
for all the nine quality indicators (see Subsection VI-B).

Fig. 4: Results on Fconvex .

Since the distribution of objective vectors in some objective
vector sets is not uniform, one may think that the approximatedobjective vector sets are far from optimal. However, as demonstrated in Subsection VI-B using Table II, some quality indicators do not prefer uniformly distributed objective vectors. For
example, as shown in Figs. 3 (j) and (h), ASLD clearly has abetter uniformity than A∆ on Fconcave . Nevertheless, Table IIshows that ∆ evaluates A∆ as better than ASLD (details are
discussed later). Such an unintuitive result may be due to ahidden property of each quality indicator. We reemphasize thata theoretical analysis of the optimal µ-distribution for m ≥ 3 isdifficult. For this reason, the (true and approximated) optimalµ-distributions for almost all quality indicators have not been
investigated on the Pareto front with m ≥ 3.
Below, we individually discuss the optimal µ-distribution
for each of the nine quality indicators (AHV , AIGD , AIGD+ ,AR2 , ANR2 , AI+ , ASE , A∆ , and APD ). Although results of
HV, IGD, and IGD+ are consistent with the previous studies,some results of other quality indicators provide insightful
information. Our findings are summarized in Section VII.
1) AHV : It should be noted that the optimal µ-distributionfor HV significantly depends on the position of the reference
vector q [4]. Objective vectors in AHV are on the edge andthe center of the Pareto front for all the front shape functions,
except for Fconvex and Fi-convex . These observations on the9(a) AHV(b) AIGD(a) AHV(b) AIGD(a) AHV(b) AIGD(c) AIGD+(d) AR2(c) AIGD+(d) AR2(c) AIGD+(d) AR2(e) ANR2(f) AI+(e) ANR2(f) AI+(e) ANR2(f) AI+(g) ASE(h) A∆(g) ASE(h) A∆(g) ASE(h) A∆(i) APD(j) ASLD(i) APD(j) ASLD(i) APD(j) ASLD

Fig. 5: Results on Fi-linear .

Fig. 6: Results on Fi-concave .

convex Pareto fronts are consistent with the previous studies[4], [12], [15]. When µ is small, AHV does not containthe extreme objective vectors even for the linear front shapefunction Flinear (e.g., AHV with µ = 10 shown in Fig. S.3 inthe supplementary file). As µ increases, AHV covers the edgeof the Pareto front on most front shape functions.
2) AIGD : The distribution of objective vectors in AIGD
is uniform on most front shape functions, even including
Fc-concave . However, objective vectors in AIGD on all front
shape functions are not on the edge of the Pareto front. Thus,
IGD may incorrectly evaluate the spread quality of objective
vector sets. The poor spread quality of AIGD is consistent
with the previous studies [3], [31].
3) AIGD+ : Our results are consistent with the resultspresented in [31]. The distribution of objective vectors in
AIGD+ is almost the same as that in AIGD on Flinear andsimilar to that in AHV on the other front shape functions. As
the µ value increases, the distributions of objective vectors in
AIGD+ and AHV become more similar. For details, see Figs.
S.3–S.50 in the supplementary file.
4) AR2 : Objective vectors in AR2 are densely distributedin the center of the Pareto front for Flinear , Fconcave , andFconvex . These results are consistent with the previous study[17] for m = 2. However, our results on Fi-linear and Fi-concave
are inconsistent with [17]. Only one objective vector in AR2
is in the center of the Pareto front for Fi-linear . Almost all

Fig. 7: Results on Fi-convex .

objective vectors in AR2 are on the edge of the Pareto front.
Similarly, all objective vectors in AR2 are on the edge of thePareto front for Fi-concave . When µ is set to a larger value,
AR2 contains a few objective vectors on the center of thePareto front for Fi-concave as shown in Figs. S.3–S.50.
Our inconsistent results on Fi-linear and Fi-concave aremainly due to the inverted triangular Pareto front. The analysis presented in [51] shows that some decomposition-based
EMOAs (e.g., MOEA/D and NSGA-III) perform poorly on
problems with inverted triangular Pareto fronts. This is because
the shape of the distribution of weight vectors used in thedecomposition-based EMOAs is different from the shape of
the Pareto front. Since R2 uses the weight vector set W thatis similar to decomposition-based EMOAs, R2 is likely to havea similar issue in decomposition-based EMOAs for m = 3.
5) ANR2 : Since NR2 was designed to approximate theHV value, ANR2 is very similar to AHV for all front shapefunctions. Although NR2 is the modified version of R2, ANR2
is dissimilar to AR2 . For example, AR2 does not contain any
objective vector in the center of the Pareto front for Fi-concave .
In contrast, ANR2 contains some objective vectors in the centerof the Pareto front for Fi-concave that is similar to AHV .
6) AI+ : As mentioned in Subsection III-B, the optimalµ-distributions for I+ on m = 2 were investigated in [19].
However, our results on m = 3 are inconsistent with [19].
The distribution of objective vectors in AI+ is not uniform

10(a) AHV(b) AIGD(c) AIGD+(d) AR2(e) ANR2(f) AI+(g) ASE(h) A∆(i) APD

Fig. 8: Approximated optimal µ-distributions on Fdisconnected .

(a) AHV(b) AIGD(c) AIGD+(d) AR2(e) ANR2(f) AI+(g) ASE(h) A∆(i) APD

Fig. 9: Approximated optimal µ-distributions on Fc-concave .
for all front shape functions (except for Fi-concave ). Thus, ourresults show that optimal µ-distributions for I+ depend on the
choice of m. As pointed out in [5], the I+ value of a given
objective vector set A depends only on one objective value of
one objective vector in A. This property of I+ may influence
its optimal µ-distributions on m = 3 in a complicated manner.
An in-depth analysis of the optimal µ-distributions for I+ on
m = 3 is another future work.
7) ASE : Similar to AIGD , the distribution of objective
vectors in ASE is uniform for all front shape functions. In
contrast, objective vectors in ASE are more widely distributed
than those in AIGD . For example, while objective vectors in
AIGD are densely distributed in the center of the Pareto front
for Fconvex , those in ASE cover the whole Pareto front. We
can conclude that SE evaluates the quality of objective vector
sets in terms of both uniformity and spread.
However, most objective vectors in ASE are on the edgeof the Pareto front (especially on Fdisconnected ). Even as theµ value increases, this characteristic can be still found (seeFigs. S.3–S.50). Thus, SE may overestimate a set of objective
vectors on the edge of the Pareto front.
8) A∆ : Since A∆ contains the extreme objective vectors
on all F , ∆ can evaluate the spread quality of objective vector
sets. In contrast, the distribution of objective vectors in A∆
is far from uniform. Thus, ∆ may incorrectly evaluate theuniformity quality of objective vector sets.
∆ uses a pairwise distance between two closest objective
vectors. As pointed out in [15], this calculation scheme in ∆can provide misleading results in some cases. This is because∆ does not consider the distance between non-closest objective
vectors. In fact, objective vectors in A∆ for Flinear seem tobe connected in a chained manner (e.g., see Fig. 2 (h)). In
contrast, SE handles the distance between all pairs of two
objective vectors. Thus, SE does not have this issue in ∆.
It should be noted that the misleading results by ∆ havealready been reported in [5], [15]. Thus, this observation is
not new. Our contribution here is that we find and analyzeA∆ on the eight front shape functions. The above-mentioned
unintuitive distribution of objective vectors in A∆ has not been
investigated in the literature.
9) APD : Unlike the other quality indicators, PD evaluates
the dissimilarity between objective vectors. For this reason,
the distribution of objective vectors in APD is not uniform on
all front shape functions. This observation is consistent withthe property of PD reported in [5], [24]. Also, APD does notcontain all the three extreme objective vectors in most cases.
Interestingly, at least two objective vectors in APD overlap.
Since some objective vectors in APD completely overlap on
Fi-concave and Fdisconnected , it looks like that µ is less than
21. Thus, PD may overestimate an objective vector set that
contains similar objective vectors. Although we calculated the
PD value using the translated version of the original sourcecode as explained in Section V, we can obtain exactly thesame results even using the original one.
To understand the influence of overlapping objective vectors
on PD, we slightly modified one of overlapping objective
vectors in APD on Flinear (Fig. 5 (i)). Two resulting objective
vector sets are denoted as Amod1
and Amod2
PD
PD . Fig. S.2 inthe supplementary file shows the distribution of objective
vectors in Amod1
and Amod2
on Flinear . The distribution ofPDPD
objective vectors in Amod1
and
Amod2
is more diverse than
PD
PD
mod1
that in APD . Nevertheless, APD and Amod2
are evaluatedPD
as worse than APD as follows: PD(APD ) = 1.74 × 105 ,mod2
5
5
PD(Amod1
PD ) = 1.32 × 10 , and PD(APD ) = 1.22 × 10 .
Recall that PD is to be maximized.
Our observation may be due to the sensitivity of PD tothe order of objective vectors. PD constructs a tree based onµ objective vectors a1 , ..., aµ in an objective vector set A.
In this procedure, each objective vector is greedily linked toits nearest unreplicated objective vector in lexical order (i.e.,from a1 to aµ ). The resulting tree depends on the order of theobjective vectors. We arranged the objective vectors in APD
on Flinear in reverse order (i.e., from aµ to a1 ). The resulting
mod3
objective vector set is denoted as Amod3
PD . Although APD
and APD contain exactly the same objective vectors, Amod3
PD

11

TABLE II: Rankings of the nine approximated optimal µ-distributions (AHV , ..., APD ) and ASLD by each quality indicator on Fconcave .
AHV

AIGD

AIGD+

AR2

ANR2

AI+

ASE

A∆

APD

ASLD

HV
IGD
IGD+
R2
NR2
I+
SE
∆
PD

1
9
2
2
2
3
7
9
7

10
1
9
10
10
5
6
3
5

3
10
1
4
3
6
4
6
8

5
7
4
1
5
4
3
4
3

2
8
3
3
1
2
5
7
9

7
2
7
8
7
1
8
5
2

4
4
5
5
4
8
1
2
6

8
6
8
7
8
9
9
1
4

9
5
10
9
9
10
10
10
1

6
3
6
6
6
7
2
8
10

Avg.

4.2

5.9

4.5

3.6

4.0

4.7

3.9

6.0

7.3

5.4

is evaluated as worse than APD in terms of PD as follows:
5
PD(Amod3
PD ) = 1.29 × 10 . Overlapping objective vectors
in APD may help to construct a larger tree by implicitly
exploiting the lexical characteristic of PD.
B. Ranking information based analysis
In Subsection VI-A, we analyzed the nine quality indicators
using the optimal µ-distributions. Here, we examine the nine
quality indicators using the ranking information.
Table II shows the rankings of the nine approximated optimal µ-distributions (AHV , ..., APD ) and ASLD by each quality
indicator on Fconcave . The average ranking values are reported
at the bottom of Table II. Tables S.2–S.6 in the supplementary
file show the rankings on Flinear , Fconvex , Fi-linear , Fi-concave ,and Fi-convex , respectively. Since ASLD cannot be generated
on Fdisconnected and Fc-concave , the rankings for Fdisconnected
and Fc-concave are not shown. Note that the quality indicator
called “SLD” does not exist. Recall that ASLD is generated
by simplex-lattice design. Due to the paper length limitation,
we mainly explain the results on Fconcave .
On the one hand, from each row in Table II, we can read theranking of the ten objective vector sets by the corresponding
quality indicator on Fconcave . For example, HV evaluates AHV
and AIGD as the best and the worst, respectively. On the other
hand, from each column in Table II, we can read the ranking of
the corresponding objective vector set by the nine indicators.
For example, ASE is ranked at fourth by HV, IGD, and NR2.
The average ranking value of each objective vector set in thebottom row of Table II (i.e., average value of each column)means how highly the corresponding objective vector set isevaluated by all the nine quality indicators. AR2 performs thebest among the 10 objective vector sets on Fconcave in terms of
the average ranking values. Thus, a good objective vector set
with respect to R2 may be highly evaluated by the other quality
indicators on Fconcave . In contrast, APD performs the worst interms of the average ranking values. This result indicates that agood objective vector set with respect to PD may be evaluatedas inferior to other objective vector sets by the other quality
indicators on Fconcave . These results are useful for selecting
a quality indicator used in indicator-based EMOAs. Note thatthe best objective vector set in terms of the average ranking
values depends on the choice of F . For example, Table S.4
shows that AR2 performs the second worst on Fi-linear .

TABLE III: Kendall τ values of the nine indicator on Fconcave .
HV
HV
1.00
IGD
-0.60
IGD+ 0.82
R2
0.78
NR2
0.96
I+
0.33
SE
0.24
∆
-0.24
PD
-0.42

IGD
-0.60
1.00
-0.69
-0.56
-0.56
-0.02
-0.02
0.20
0.20

IGD+
0.82
-0.69
1.00
0.69
0.78
0.33
0.33
-0.16
-0.42

R2
0.78
-0.56
0.69
1.00
0.73
0.29
0.29
-0.11
-0.29

NR2
0.96
-0.56
0.78
0.73
1.00
0.38
0.29
-0.20
-0.47

I+
0.33
-0.02
0.33
0.29
0.38
1.00
-0.07
-0.11
-0.11

SE
0.24
-0.02
0.33
0.29
0.29
-0.07
1.00
0.16
-0.47

∆
-0.24
0.20
-0.16
-0.11
-0.20
-0.11
0.16
1.00
0.11

PD
-0.42
0.20
-0.42
-0.29
-0.47
-0.11
-0.47
0.11
1.00

Each quality indicator I evaluates its optimal µ-distribution
AI as the best. For this reason, all diagonal elements in Table
II are “1”. The same results can be found in the rankings
on the other front shape functions, except for the ranking of
AIGD+ by IGD+ on Flinear in Table S.2 and ANR2 by NR2
on Fi-concave in Table S.5. IGD+ and NR2 evaluate AIGD and
AHV as the best in Tables S.2 and S.5, respectively. These
exceptional results are due to the similarity between IGD and
IGD+ , and the similarity between HV and NR2.
As in [38], Table III shows the rank-based nonlinear Kendall
rank correlation τ values of the nine quality indicator on
Fconcave . Tables S.7–S.11 in the supplementary file show theτ values of the nine quality indicator on Flinear , Fconvex ,Fi-linear , Fi-concave , and Fi-convex , respectively. The τ valuein Table III represents the similarity of the rankings by two
quality indicators I1 and I2 on the 10 objective vector sets
(AHV , ..., ASLD ). The τ values in Table III are symmetric.
The range of the τ value is [−1, 1]. While the non-negative
τ value means that I1 and I2 are consistent with each other,
the negative τ value means that I1 and I2 are inconsistent
with each other. For example, Table III shows that HV isinconsistent with IGD (τ = −0.60), but HV is consistent with
IGD+ (τ = 0.82).
As seen from Table II, no optimal µ-distribution is ranked
as the best by multiple quality indicators due to the conflicting
properties of the nine quality indicators. For example, AIGD
is ranked as the best by IGD but the worst by HV, R2, andNR2. In contrast, the rankings by HV and NR2 are almost thesame on all front shape functions (except for the rankings of
AHV and ANR2 ). In fact, the τ value of HV and NR2 in Table
III is 0.96, which indicates HV and NR2 are highly correlatedwith each other. Although the rankings by IGD and IGD+

12

are almost the same on the linear Pareto fronts (Flinear andFi-linear ), they are dissimilar on the non-linear Pareto fronts(Fconcave , Fconvex , Fi-concave , and Fi-convex ). The τ value of
IGD and IGD+ in Table III is also −0.60 on Fconcave . These
results of IGD and IGD+ are consistent with [31], [40].
The simplex lattice-based ASLD can be viewed as the
optimal objective vector set obtained by decomposition-based
EMOAs. Since ASLD has a good uniformity, ASLD is evaluated as the second best by SE on Fconcave . In contrast,
ASLD is evaluated as “poor” by the other quality indicators
in most cases (except for the results on Fi-linear in Table
S.4). In fact, the distribution of objective vectors in ASLD
is not identical to optimal µ-distributions for the nine quality
indicators (see Figs. 2–7). The poor results of ASLD suggestthat decomposition-based EMOAs need an adaptive weightvector method considering the optimal µ-distributions for agiven quality indicator.
C. Analysis for m ≥ 4
Since it is almost impossible to clearly show the distributionof objective vectors of a problem with m ≥ 4 in an understandable manner [56], [57], this paper focused on m = 3. We
reemphasize that our analysis for m = 3 itself is already an
important contribution to the EMO community. Nevertheless,
it is interesting to investigate whether the results for m = 3
can be generalized to the case of m ≥ 4.
Since the distribution of objective vectors for m ≥ 4 isunclear, our analysis here is based on the ranking information
as in Subsection VI-B. Tables S.12–S.23 in the supplementary
file show the rankings of the nine approximated optimal µdistributions (AHV , ..., APD ) with µ = 21 by each quality indicator on six front shape functions (Flinear , Fconcave ,Fconvex , Fi-linear , Fi-concave , and Fi-convex ) with m = 5 and8, respectively. Since ASLD with µ = 21 cannot be generated
for m = 5 and 8, ASLD is removed from this experiment. The
size of the reference point set R was set to 3 060 for m = 5
and 5 148 for m = 8. We generated R for m = 8 using
the two-layered version of simplex-lattice design [49]. Tables
S.24–S.35 in the supplementary file also show the Kendall
rank correlation τ values of the nine quality indicators on thesix front shape functions with m = 5 and 8, respectively.
Although most results for m ≥ 4 are consistent with theresults for m = 3, some exceptional results can be found.
ASE is ranked as the best by ∆ on Flinear and Fi-linearwith m = 8, and Fconvex and Fi-concave with m ∈ {5, 8}.
This means that a better approximated optimal µ-distributionfor ∆ can be found by optimizing SE. Similarly, AR2 isranked as the best by I+ on Fi-concave with m ∈ {5, 8}.
This may be because of the problem difficulty in finding
the optimal µ-distribution for ∆ and I+ . For example, as
mentioned in Subsection VI, the I+ value of A depends only
on one objective value of one objective vector in A. As in
the Schwefel 2.21 function f (x) = maxi∈{1,...,n} {|xi |} [58],this property of I+ can incorporate a strong nonseparability
between variables θ1 , ..., θd in the proposed formulation. In
addition, this property can make some areas of the fitness
landscape of the proposed formulation plateaus. These problem difficulties become pronounced with the increase of m.

TABLE IV: Best objective vector sets on each front shape function
with m ∈ {3, 5, 8} in terms of the average ranking.
Flinear

Fconcave

FconvexFi-linearFi-concave

Fi-convex

3

AHV , ANR2

AR2

AHV

AHV

AHV

AHV

5

ANR2 , A∆

ANR2

AHV

ANR2

AHV

AIGD+

8

A∆

AIGD+

AIGD

AIGD+

AHV

AIGD+

A landscape analysis of the proposed formulation for eachquality indicator is needed for deeper understanding.
Table IV shows the best objective vector sets on each F withm ∈ {3, 5, 8} in terms of the average ranking. We can see thatthe best objective vector set in terms of the average ranking
depends on m in addition to the choice of F . For example,
AR2 , ANR2 , and AIGD+ perform the best on Fconcave with
m = 3, 5, and 8, respectively. These results indicate thenecessity of an adaptive indicator selection mechanism inindicator-based EMOAs (e.g., [34]).
In addition, the correlation between HV and NR2 becomes
weak with the increase of m. For example, the τ values of HV
and NR2 on Flinear are 0.96 for m = 3, 0.89 for m = 5, and0.22 for m = 8 (see Tables III, S.24, and S.30, respectively).
This observation indicates that the approximation performanceof NR2 is likely to deteriorate with the increase of m.
In summary, the results here show that some results form = 3 cannot be always generalized to the case of m ≥ 4.
Thus, our observations are only for m = 3 as emphasizedin the title of this paper. Further analysis for many-objective
problems is need in future research.
D. Analysis of a unary version of an M -nary indicator
As in I+ , some binary or M -nary quality indicators canbe converted into their unary versions by using the reference
vector set R as the compared vector set. However, it is questionable that such a converted unary quality indicator works
well. Here, we demonstrate that the proposed formulation canbe useful to analyze the unary version of an M -nary indicator.
Since a number of M -nary quality indicators have been
proposed, it is very difficult to investigate all of them. Instead,
we selected the diversity comparison indicator (DCI) [59]since it is recommended in [5]. DCI relatively compares M
objective vector sets A1 , ..., AM . A large DCI value indicatesthat the corresponding A has a good diversity among A1 , ...,
AM . We explain DCI in Section S.4 in the supplementary file.
For the DCI calculation, we used the source code provided by
the authors of [59].
We converted DCI into the unary indicator using the reference vector set R. To generate R, we used the same method
described in Section V. Thus, the DCI value of A is obtainedby comparing A to R by DCI. As recommend in [59], we set
div in DCI to 19 for m = 3 in this study. Since we want toexamine the general property of DCI, we did not finely tune
div for each F . We set µ to 21. We investigate the influence ofthe size of R on DCI. Below, we denote DCI with |R| = 21
and |R| = 28 as DCI21 and DCI28, respectively.

13(a) ADCI21(b) ADCI28

Fig. 10: ADCI21 and ADCI28 on Flinear .

Fig. 10 shows the distribution of objective vectors in
ADCI21 and ADCI28 on Flinear . Figs. S.51–S.55 in the supplementary file show ADCI21 (and ADCI28 ) on Fconcave ,Fconvex , Fi-linear , Fi-concave , and Fi-convex , respectively. On
the one hand, Fig. 10 (a) shows that the objective vectors in ADCI21 are uniformly distributed. Thus, DCI21 is
likely to prefer uniformly distributed objective vectors. In
fact, DCI21(ADCI21 ) = 1 and DCI21(ASLD ) = 1, where
ASLD is shown in Fig. 2 (j). On the other hand, Fig. 10(b) shows that the distribution of the objective vectors in
ADCI28 is not uniform. This result means that DCI28 does
not prefer uniformly distributed objective vectors. Actually,
DCI28 evaluates ADCI28 as better than ASLD as follows:
DCI28(ADCI28 ) = 0.75 and DCI28(ASLD ) = 0.67.
In summary, we can conclude that it is possible to design aunary version of DCI, but its result is sensitive to the size of R.
In the unary version of DCI, the size of R and µ should be thesame. Otherwise, non-uniformly distributed objective vectors
are highly evaluated by DCI (as shown in Fig. 10 (b)).
VII. C ONCLUSION
We have analyzed the nine quality indicators using theirapproximated optimal µ-distributions for m = 3. First, we
proposed the problem formulation of finding the optimal µdistribution for each quality indicator on Pareto fronts for arbitrary number of objectives (m ≥ 2). Then, we approximatedthe optimal µ-distributions for the nine quality indicators on
the Pareto fronts of eight problems with m = 3. We analyzed
the nine quality indicators based on the optimal µ-distribution
based approach and the ranking information based approach.
We also examined the generality of results for m = 3 withrespect to m and a unary version of an M -nary indicator.
Our observations for m = 3 can be summarized as follows:i) Objective vectors in the approximated optimal µdistributions for R2 are biased on the edge of the Pareto
front for Fi-linear and Fi-concave . These results for m = 3
are inconsistent with the results for m = 2.
ii) The approximated optimal µ-distributions for NR2 andHV are similar. NR2 and HV are correlated with eachother in terms of the Kendall rank correlation (τ = 0.96
on Fconcave ). Thus, NR2 can substitute for HV.
iii) While the optimal µ-distributions for I+ for m = 2 arealmost uniform, those for m = 3 are not uniform evenon the linear front.
iv) Although SE can evaluate both uniformity and spreadquality of objective vector sets, SE prefers objective
vectors that are on the edge of the Pareto front.

v) Since the optimal µ-distributions for ∆ are far from
uniform, ∆ may incorrectly evaluate the uniformity
quality of objective vector sets.
vi) PD may overestimate an objective vector set with asmall dissimilarity. PD is also sensitive to the order ofobjective vectors.
vii) Most quality indicators evaluate objective vectors
generated by simplex-lattice design as poor. Thus,
decomposition-based EMOAs need an adaptive weightvector method when they try to search for a good objective vector set with respect to each quality indicator.
viii) A unary version of DCI is sensitive to the size of thereference vector set.
We emphasize that the above-mentioned observations on
the Pareto front for m = 3 could not be obtained without
searching for the optimal distributions of objective vectors by
the proposed problem formulation. As shown in Subsection
VI-C, our results for m = 3 cannot be always generalized
to the case of m ≥ 4. Further analysis is needed for manyobjective problems in future research.
A further analysis of the optimal µ-distributions for R2 andI+ on the Pareto front for m = 3 is needed. The undesirable
property of PD could be addressed by using an algorithm forfinding a minimum spanning tree (e.g., [60]). The problemof finding the optimal µ-distribution in the proposed problem
formulation itself can be viewed as a “real-world” black-box
single-objective optimization problem. Although we used LSHADE in this paper, we believe that benchmarking various
single-objective optimizers on the proposed problem formulation is beneficial for both the EMO community and theevolutionary single-objective optimization community.
ACKNOWLEDGMENT
This work was supported by the Program for Guangdong Introducing Innovative and Enterpreneurial Teams(Grant No. 2017ZT07X386), Shenzhen Peacock Plan (Grant
No. KQTD2016112514355531), the Science and Technology Innovation Committee Foundation of Shenzhen (Grant
No. ZDSYS201703031748284), the Program for University Key Laboratory of Guangdong Province (Grant No.
2017KSYS008), and National Natural Science Foundation ofChina (Grant No. 61876075).
R EFERENCES[1] K. Miettinen, Nonlinear Multiobjective Optimization. Springer, 1998.
[2] K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms.
John Wiley & Sons, 2001.
[3] H. Ishibuchi, R. Imada, Y. Setoguchi, and Y. Nojima, “Reference Point
Specification in Inverted Generational Distance for Triangular Linear
Pareto Front,” IEEE TEVC, vol. 22, no. 6, pp. 961–975, 2018.
[4] ——, “How to Specify a Reference Point in Hypervolume Calculation
for Fair Performance Comparison,” Evol. Comput., vol. 26, no. 3, 2018.
[5] M. Li and X. Yao, “Quality Evaluation of Solution Sets in Multiobjective
Optimisation: A Survey,” ACM CSUR, vol. 52, no. 2, pp. 26:1–26:38,
2019.
[6] J. Knowles and D. Corne, “On Metrics for Comparing NondominatedSets,” in IEEE CEC, 2002, pp. 711–716.
[7] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and V. G. da Fonseca, “Performance assessment of multiobjective optimizers: an analysis
and review,” IEEE TEVC, vol. 7, no. 2, pp. 117–132, 2003.
[8] T. Okabe, Y. Jin, and B. Sendhoff, “A Critical Survey of Performance

14

Indices for Multi-Objective Optimisation,” in IEEE CEC, 2002, pp. 878–
885.
[9] E. Zitzler and L. Thiele, “Multiobjective Optimization Using Evolutionary Algorithms - A Comparative Case Study,” in PPSN, 1998, pp.
292–304.
[10] C. A. C. Coello and M. R. Sierra, “A Study of the Parallelization ofa Coevolutionary Multi-objective Evolutionary Algorithm,” in MICAI,
2004, pp. 688–697.
[11] H. Ishibuchi and T. Murata, “A multi-objective genetic local searchalgorithm and its application to flowshop scheduling,” IEEE Trans. SMC,
Part C, vol. 28, no. 3, pp. 392–403, 1998.
[12] A. Auger, J. Bader, D. Brockhoff, and E. Zitzler, “Theory of thehypervolume indicator: optimal µ-distributions and the choice of thereference point,” in FOGA, 2009, pp. 87–102.
[13] K. Deb, S. Agrawal, A. Pratap, and T. Meyarivan, “A fast and elitist
multiobjective genetic algorithm: NSGA-II,” IEEE TEVC, vol. 6, no. 2,
pp. 182–197, 2002.
[14] Y. Wang, L. Wu, and X. Yuan, “Multi-objective self-adaptive differential
evolution with elitist archive and crowding entropy-based diversity
measure,” Soft Comput., vol. 14, no. 3, pp. 193–209, 2010.
[15] S. Jiang, Y. Ong, J. Zhang, and L. Feng, “Consistencies and Contradictions of Performance Metrics in Multiobjective Optimization,” IEEE
Trans. Cyber., vol. 44, no. 12, pp. 2391–2404, 2014.
[16] A. Auger, J. Bader, D. Brockhoff, and E. Zitzler, “Hypervolumebased multiobjective optimization: Theoretical foundations and practical
implications,” Theor. Comput. Sci., vol. 425, pp. 75–103, 2012.
[17] D. Brockhoff, T. Wagner, and H. Trautmann, “On the properties of theR2 indicator,” in GECCO, 2012, pp. 465–472.
[18] T. Glasmachers, “Optimized Approximation Sets for Low-Dimensional
Benchmark Pareto Fronts,” in PPSN, 2014, pp. 569–578.
[19] K. Bringmann, T. Friedrich, and P. Klitzke, “Efficient computationof two-dimensional solution sets maximizing the epsilon-indicator,” inIEEE CEC, 2015, pp. 970–977.
[20] N. Hansen and A. Ostermeier, “Completely derandomized selfadaptation in evolution strategies,” Evol. Comput., vol. 9, no. 2, pp.
159–195, 2001.
[21] H. Ishibuchi, H. Masuda, Y. Tanigaki, and Y. Nojima, “Modified
Distance Calculation in Generational Distance and Inverted Generational
Distance,” in EMO, 2015, pp. 110–125.
[22] K. Shang, H. Ishibuchi, M. Zhang, and Y. Liu, “A new R2 indicator for
better hypervolume approximation,” in GECCO, 2018, pp. 745–752.
[23] D. P. Hardin and E. B. Saff, “Discretizing Manifolds via Minimum
Energy Points,” Notices of the AMS, vol. 51, no. 10, pp. 1186–1194,
2004.
[24] H. Wang, Y. Jin, and X. Yao, “Diversity Assessment in Many-Objective
Optimization,” IEEE Trans. Cyber., vol. 47, no. 6, pp. 1510–1522, 2017.
[25] E. Zitzler, L. Thiele, and J. Bader, “On Set-Based Multiobjective
Optimization,” IEEE TEVC, vol. 14, no. 1, pp. 58–79, 2010.
[26] J. Knowles, L. Thiele, and E. Zitzler, “A tutorial on the performanceassessment of stochastic multiobjective optimizers,” ETH Zurich, Tech.
Rep. TIK-214, 2006.
[27] D. A. V. Veldhuizen and G. B. Lamont, “Multiobjective Evolutionary
Algorithm Research: A History and Analysis,” Air Force Institute ofTechnology, Tech. Rep. TR-98-03, 1998.
[28] D. A. V. Veldhuizen, “Multiobjective evolutionary algorithms: Classifications, analyses, and new innovations,” Air University, Tech. Rep.
AFIT/DS/ENG/99-01, 1999.
[29] E. Zitzler, D. Brockhoff, and L. Thiele, “The Hypervolume Indicator
Revisited: On the Design of Pareto-compliant Indicators Via Weighted
Integration,” in EMO, 2006, pp. 862–876.
[30] O. Schütze, X. Esquivel, A. Lara, and C. A. C. Coello, “Using the
Averaged Hausdorff Distance as a Performance Measure in Evolutionary
Multiobjective Optimization,” IEEE TEVC, vol. 16, no. 4, pp. 504–522,
2012.
[31] H. Ishibuchi, R. Imada, N. Masuyama, and Y. Nojima, “Comparisonof Hypervolume, IGD and IGD+ from the Viewpoint of OptimalDistributions of Solutions,” in EMO, 2019, pp. 332–345.
[32] M. P. Hansen and A. Jaszkiewicz, “Evaluating the quality of approximations to the non-dominated set,” Poznan University of Technology,
Tech. Rep. IMM-REP-1998-7, 1998.
[33] H. Ishibuchi, N. Tsukamoto, Y. Sakane, and Y. Nojima, “Indicator-based
evolutionary algorithm with hypervolume approximation by achievement
scalarizing functions,” in GECCO, 2010, pp. 527–534.
[34] J. G. Falcón-Cardona and C. A. C. Coello, “A multi-objective evolutionary hyper-heuristic based on multiple indicator-based density estimators,” in GECCO, 2018, pp. 633–640.
[35] A. R. Solow and S. Polasky, “Measuring biological diversity,” Environ.

Ecol. Stat., vol. 1, no. 2, pp. 95–103, 1994.
[36] M. Ravber, M. Mernik, and M. Crepinsek, “The impact of Quality
Indicators on the rating of Multi-objective Evolutionary Algorithms,”
Appl. Soft Comput., vol. 55, pp. 265–275, 2017.
[37] S. Wessing and B. Naujoks, “Sequential parameter optimization for
multi-objective problems,” in IEEE CEC, 2010, pp. 1–8.
[38] A. Liefooghe and B. Derbel, “A Correlation Analysis of Set Quality
Indicator Values in Multiobjective Optimization,” in GECCO, 2016, pp.
581–588.
[39] H. Ishibuchi, H. Masuda, and Y. Nojima, “A Study on Performance Evaluation Ability of a Modified Inverted Generational Distance Indicator,”
in GECCO, 2015, pp. 695–702.
[40] L. C. T. Bezerra, M. López-Ibáñez, and T. Stützle, “An Empirical
Assessment of the Properties of Inverted Generational Distance on
Multi- and Many-Objective Optimization,” in EMO, 2017, pp. 31–45.
[41] K. Deb, L. Thiele, M. Laumanns, and E. Zitzler, “Scalable Test Problems for Evolutionary Multi-Objective Optimization,” in Evolutionary
Multiobjective Optimization. Theoretical Advances and Applications.
Springer, 2005, pp. 105–145.
[42] E. Zitzler, K. Deb, and L. Thiele, “Comparison of Multiobjective
Evolutionary Algorithms: Empirical Results,” Evol. Comput., vol. 8,
no. 2, pp. 173–195, 2000.
[43] K. Bringmann, T. Friedrich, and P. Klitzke, “Two-dimensional subset
selection for hypervolume and epsilon-indicator,” in GECCO, 2014, pp.
589–596.
[44] N. Beume, B. Naujoks, and M. T. M. Emmerich, “SMS-EMOA:
multiobjective selection based on dominated hypervolume,” EJOR, vol.
181, no. 3, pp. 1653–1669, 2007.
[45] S. Z. Martı́nez, C. A. C. Coello, H. E. Aguirre, and K. Tanaka, “A
Review of Features and Limitations of Existing Scalable Multiobjective
Test Suites,” IEEE TEVC, vol. 23, no. 1, pp. 130–142, 2019.
[46] A. Jaszkiewicz, “On the performance of multiple-objective genetic local
search on the 0/1 knapsack problem - a comparative experiment,” IEEE
TEVC, vol. 6, no. 4, pp. 402–412, 2002.
[47] H. Ishibuchi, H. Masuda, and Y. Nojima, “Pareto Fronts of ManyObjective Degenerate Test Problems,” IEEE TEVC, vol. 20, no. 5, pp.
807–813, 2016.
[48] Y. Tian, X. Xiang, X. Zhang, R. Cheng, and Y. Jin, “Sampling Reference
Points on the Pareto Fronts of Benchmark Multi-Objective Optimization
Problems,” in IEEE CEC, 2018, pp. 1–6.
[49] K. Deb and H. Jain, “An evolutionary many-objective optimization
algorithm using reference-point-based nondominated sorting approach,part I: solving problems with box constraints,” IEEE TEVC, vol. 18,
no. 4, pp. 577–601, 2014.
[50] H. Jain and K. 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 TEVC, vol. 18, no. 4, pp. 602–622, 2014.
[51] H. Ishibuchi, Y. Setoguchi, H. Masuda, and Y. Nojima, “Performanceof Decomposition-Based Many-Objective Algorithms Strongly Depends
on Pareto Front Shapes,” IEEE TEVC, vol. 21, no. 2, pp. 169–190, 2017.
[52] R. Tanabe and A. S. Fukunaga, “Improving the search performance of
SHADE using linear population size reduction,” in IEEE CEC, 2014,
pp. 1658–1665.
[53] R. Storn and K. Price, “Differential Evolution - A Simple and Efficient
Heuristic for Global Optimization over Continuous Spaces,” J. Global
Optimiz., vol. 11, no. 4, pp. 341–359, 1997.
[54] R. L. While, L. Bradstreet, and L. Barone, “A Fast Way of Calculating
Exact Hypervolumes,” IEEE TEVC, vol. 16, no. 1, pp. 86–95, 2012.
[55] I. Das and J. E. Dennis, “Normal-Boundary Intersection: A New
Method for Generating the Pareto Surface in Nonlinear Multicriteria
Optimization Problems,” SIAM J. Optimiz, vol. 8, no. 3, pp. 631–657,
1998.
[56] T. Tusar and B. Filipic, “Visualization of Pareto Front Approximations
in Evolutionary Multiobjective Optimization: A Critical Review and theProsection Method,” IEEE TEVC, vol. 19, no. 2, pp. 225–245, 2015.
[57] M. Li, L. Zhen, and X. Yao, “How to Read Many-Objective Solution
Sets in Parallel Coordinates [Educational Forum],” IEEE CIM, vol. 12,
no. 4, pp. 88–100, 2017.
[58] X. Yao, Y. Liu, and G. Lin, “Evolutionary Programming Made Faster,”IEEE TEVC, vol. 3, no. 2, pp. 82–102, 1999.
[59] M. Li, S. Yang, and X. Liu, “Diversity Comparison of Pareto Front
Approximations in Many-Objective Optimization,” IEEE Trans. Cyber.,vol. 44, no. 12, pp. 2568–2584, 2014.
[60] B. Chazelle, “A minimum spanning tree algorithm with InverseAckermann type complexity,” J. ACM, vol. 47, no. 6, pp. 1028–1047,
2000.

15

Ryoji Tanabe is currently an Assistant Professor
with Faculty of Environment and Information Sciences, Yokohama National University, Japan. He
was a Research Assistant Professor with Departmentof Computer Science and Engineering, Southern
University of Science and Technology, China, from2017 to 2019. He was a Post-Doctoral Researcher
with ISAS/JAXA, Japan, from 2016 to 2017. Hereceived his Ph.D. in Science from The University ofTokyo, Japan, in 2016. His research interests include
stochastic single- and multi-objective optimization
algorithms, parameter control in evolutionary algorithms, and automaticalgorithm configuration.

Hisao Ishibuchi (M’93-SM’10-F’14) received the
B.S. and M.S. degrees in precision mechanics from
Kyoto University, Kyoto, Japan, in 1985 and 1987,
respectively, and the Ph.D. degree in computer science from Osaka Prefecture University, Sakai, Osaka, Japan, in 1992. He was with Osaka Prefecture University in 1987-2017. Since 2017, he isa Chair Professor at Southern University of Science and Technology, China. His research interests include fuzzy rule-based classifiers, evolutionary multi-objective and many-objective optimization,
memetic algorithms, and evolutionary games.
Dr. Ishibuchi was the IEEE Computational Intelligence Society (CIS) VicePresident for Technical Activities (2010-2013), an AdCom member of theIEEE CIS (2014-2019), and the Editor-in-Chief of the IEEE COMPUTATIONAL INTELLIGENCE MAGAZINE (2014-2019).