Bezier Simplex Fitting: Describing Pareto Fronts of Simplicial Problems with Small Samples in Multi-objective Optimization

TL;DR

This paper introduces a Bezier simplex model for efficiently approximating Pareto fronts in multi-objective optimization, reducing sample size requirements by exploiting the geometric structure of solution sets.

Contribution

The paper proposes a novel Bezier simplex fitting method that decomposes high-dimensional surface fitting into low-dimensional tasks, with proven approximation capabilities.

Findings

Accurately approximates high-dimensional Pareto fronts with small samples

Demonstrates effectiveness on synthetic and real-world problems

Enables better trade-off analysis in post-optimization

Abstract

Multi-objective optimization problems require simultaneously optimizing two or more objective functions. Many studies have reported that the solution set of an M-objective optimization problem often forms an (M-1)-dimensional topological simplex (a curved line for M=2, a curved triangle for M=3, a curved tetrahedron for M=4, etc.). Since the dimensionality of the solution set increases as the number of objectives grows, an exponentially large sample size is needed to cover the solution set. To reduce the required sample size, this paper proposes a Bezier simplex model and its fitting algorithm. These techniques can exploit the simplex structure of the solution set and decompose a high-dimensional surface fitting task into a sequence of low-dimensional ones. An approximation theorem of Bezier simplices is proven. Numerical experiments with synthetic and real-world optimization problems demonstrate that the proposed method achieves an accurate approximation of high-dimensional solution sets with small samples. In practice, such an approximation will be conducted in the post-optimization process and enable a better trade-off analysis.

Tables

Table 1. Table 1: GD and IGD (avg.  ± plus-or-minus \pm s.d. over 20 trials) with subsample size ( 1 , 3 ) 1 3 (1,3) for Schaffer, ConstrEx and Osyczka2; ( 1 , 2 , 1 ) 1 2 1 (1,2,1) for 3-MED, Viennet2 and 5-MED; ( 1 , 5 , 5 ) 1 5 5 (1,5,5) for S3TD. The best scores are shown in bold. ( : ∗ ∗ p < 0.05 , : ∗ p < 0.1 {}^{**}:p<0.05,{}^{*}:p<0.1 )

Problem	Inductive skeleton			All-at-once			Response surface
	Iteration	GD	IGD	Iteration	GD	IGD	GD	IGD
Schaffer	3.00 $\pm$ 0.00	2.50e-10 $\pm$ 3.76e-11	2.49e-02 $\pm$ 2.49e-12	3.00 $\pm$ 0.00	2.51e-10 $\pm$ 3.71e-11	2.49e-02 $\pm$ 2.52e-12	5.26e-01 $\pm$ 1.78e+00	8.13e-02 $\pm$ 5.54e-02
ConstrEx	3.00 $\pm$ 0.00	2.33e-02 $\pm$ 2.32e-02	4.14e-02 $\pm$ 1.24e-02	3.00 $\pm$ 0.00	2.34e-02 $\pm$ 2.32e-02	4.13e-02 $\pm$ 1.24e-02	1.47e-02 $\pm$ 9.15e-03∗∗	2.48e-02 $\pm$ 9.24e-03∗∗
Osyczka2	3.45 $\pm$ 1.12	6.02e-02 $\pm$ 2.72e-02	8.27e-02 $\pm$ 3.49e-02	3.20 $\pm$ 0.51	6.08e-02 $\pm$ 2.86e-02	8.33e-02 $\pm$ 3.70e-02	2.79e-01 $\pm$ 3.34e-01	1.01e-01 $\pm$ 4.22e-02
3-MED	2.45 $\pm$ 0.50	3.99e-01 $\pm$ 9.77e-01∗∗	6.16e-02 $\pm$ 2.51e-02∗∗	3.00 $\pm$ 0.00	1.02e+00 $\pm$ 2.84e+00	1.05e-01 $\pm$ 6.43e-02	1.39e+00 $\pm$ 2.38e-01	6.75e-02 $\pm$ 5.37e-03
Viennet2	2.65 $\pm$ 0.48	2.51e+00 $\pm$ 3.73e+00∗∗	6.47e-02 $\pm$ 9.91e-03	3.68 $\pm$ 0.98	2.24e+09 $\pm$ 9.40e+09	2.42e-01 $\pm$ 1.02e-01	3.10e+00 $\pm$ 2.58e+00	6.89e-02 $\pm$ 1.73e-02
5-MED	2.60 $\pm$ 0.49	2.55e-01 $\pm$ 3.13e-01∗∗	7.94e-02 $\pm$ 1.40e-02∗∗	3.00 $\pm$ 0.00	4.19e+00 $\pm$ 1.24e+01	1.66e-01 $\pm$ 4.01e-02	4.89e+00 $\pm$ 7.52e-01	1.02e-01 $\pm$ 4.89e-03
S3TD	3.25 $\pm$ 1.58	2.04e-01 $\pm$ 4.04e-02∗	1.37e-01 $\pm$ 3.15e-02	36.20 $\pm$ 33.81	5.58e-01 $\pm$ 5.99e-01	6.31e-01 $\pm$ 6.36e-01	7.17e-01 $\pm$ 2.14e-01	9.88e-02 $\pm$ 6.67e-03∗∗

Table 2. Table 2: GD and IGD (avg.  ± plus-or-minus \pm s.d. over 20 trials) with sample size ( 1 , 2 , 1 ) 1 2 1 (1,2,1) for the graph of 5-MED.

Inductive skeleton			All-at-once
Iteration	GD	IGD	Iteration	GD	IGD
2.84 $\pm$ 0.36	3.36e-1 $\pm$ 9.48e-2	2.05e-1 $\pm$ 3.45e-2	3.00 $\pm$ 0.00	1.38e-00 $\pm$ 2.00e-00	2.79e-01 $\pm$ 8.21e-02

Table 3. Table 3: Sample size on the Pareto front for each subproblem.

Problem	1-obj	2-obj	3-obj
Schaffer	1 $(\times 2)$	201 $(\times 1)$	–
ConstrEx	1 $(\times 2)$	679 $(\times 1)$	–
Osyczka2	1 $(\times 2)$	570 $(\times 1)$	–
3-MED	1 $(\times 3)$	17 $(\times 3)$	153 $(\times 1)$
Viennet2	1 $(\times 3)$	811 $(f_1,f_2)$, 1310 $(f_2,f_3)$, 1618 $(f_3,f_1)$	8122 $(\times 1)$
5-MED	1 $(\times 5)$	15 $(\times 10)$	105 $(\times 10)$
S3TD	1 $(\times 4)$	19 $(f_1,f_2)$, 9 $(f_1,f_3)$, 21 $(f_1,f_4)$, 6 $(f_2,f_3)$, 13 $(f_2,f_4)$, 25 $(f_3,f_4)$	37 $(f_1,f_2,f_3)$, 24 $(f_1,f_2,f_4)$, 25 $(f_1,f_3,f_4)$, 30 $(f_2,f_3,f_4)$

Equations

$$\mathbf{f}(\mathbf{x}):=(f_1(\mathbf{x}),\dots ,f_M(\mathbf{x}))$$

$$\mathbf{x}\in X(\subseteq {\mathbb{R}}^L).$$

$$\displaystyle\bm{f}(\bm{x})\prec\bm{f}(\bm{y})\xLeftrightarrow$$

$$f_m(\mathbf{x})<f_m(\mathbf{y})\text{for some}m=1,\dots ,M,$$

$$X^{\ast }(\mathbf{f}):=\left\{\mathbf{x}\in X\right\}f(\mathbf{y})\nprec f(\mathbf{x})\text{for all}\mathbf{y}\in X$$

$$\mathbf{f}{X}^{\ast }(\mathbf{f}):=\left\{\mathbf{f}(\mathbf{x})\in {\mathbb{R}}^M\right\}\mathbf{x}\in X^{\ast }(\mathbf{f})$$

$${\Delta }^{M-1}:=\left\{(t_1,\dots ,t_M)\in {\mathbb{R}}^M\right\}0\le t_m,\sum _{m\in M}{t}_m=1$$

$${\Delta }^J:=\left\{(t_1,\dots ,t_M)\in {\Delta }^{M-1}\right\}{t}_m=0(m\notin J)$$

$${\mathbf{f}}_J:=(f_i{)}_{i\in J}:X\to {\mathbb{R}}^{|J|}$$

$${\Delta }^{(m)}:=\bigcup _{J\subseteq I,|J|=m}{\Delta }^J$$

$$\mathbf{\varphi }{|}_{{\Delta }^J}:{\Delta }^J$$

$$\mathbf{f}\circ \mathbf{\varphi }{|}_{{\Delta }^J}:{\Delta }^J$$

$$X^{\ast (m)}$$

$$\mathbf{f}{X}^{\ast (m)}$$

$$\mathbf{b}(t):=\sum _{d=0}^D\left(\genfrac{}{}{0ex}{}{D}{d}\right)t^d(1-t{)}^{(D-d)}{\mathbf{p}}_d(0\le t\le 1),$$

$$\begin{array}{rl}\underset{t_n,{\mathbf{p}}_d}{\mathrm{minimize}}& \sum _{n=1}^N{\Vert \mathbf{b}(t_n)-{\mathbf{x}}_n\Vert }^2\\ \text{subject to }& 0\le t_n\le 1(n=1,\dots ,N)\end{array}$$

$$\left<\frac{\partial}{\partial t}\Big{|}_{t=t_{n}}\bm{b}(t),\ \bm{b}(t_{n})-\bm{x}_{n}\right>=0.$$

$${\mathbb{N}}_D^M:=\left\{(d_1,\dots ,d_M)\in {\mathbb{N}}^M\right\}\sum _{m=1}^M{d}_m=D.$$

$$\mathbf{b}(\mathbf{t}):=\sum _{\mathbf{d}\in {\mathbb{N}}_D^M}\left(\genfrac{}{}{0ex}{}{D}{\mathbf{d}}\right){\mathbf{t}}^{\mathbf{d}}{\mathbf{p}}_{\mathbf{d}},$$

$$\left(\genfrac{}{}{0ex}{}{D}{\mathbf{d}}\right):=\frac{D!}{d_1!d_2!\cdots d_M!}.$$

$$\underset{i\to \infty }{\lim }\underset{\mathbf{t}\in {\Delta }^{M-1}}{\sup }|\mathbf{\varphi }(\mathbf{t})-{\mathbf{b}}^{(i)}(\mathbf{t})|=0.$$

$$G^{\ast }(\mathbf{f})=\left\{(\mathbf{x},\mathbf{f}(x))\in {\mathbb{R}}^L\times {\mathbb{R}}^M\right\}\mathbf{x}\in X^{\ast }(\mathbf{f}),$$

$$\mathbf{\varphi }\times (\mathbf{f}\circ \mathbf{\varphi }):{\Delta }^{M-1}\to G^{\ast }(\mathbf{f})$$

$$\underset{i\to \infty }{\lim }{d}_{\mathrm{H}}(X^{\ast },B^{(i)})=0$$

$$\begin{array}{rl}\underset{{\mathbf{t}}_n,{\mathbf{p}}_{\mathbf{d}}}{\mathrm{minimize}}& \sum _{n=1}^N{\Vert \mathbf{b}({\mathbf{t}}_n)-{\mathbf{x}}_n\Vert }^2\\ \text{subject to }& {\mathbf{t}}_n\in {\Delta }^{M-1}(n=1,\dots ,N)\end{array}$$

$$\left<\frac{\partial}{\partial t_{m}}\Big{|}_{\bm{t}=\bm{t}_{n}}\bm{b}(\bm{t}),\ \bm{b}(\bm{t}_{n})-\bm{x}_{n}\right>=0~{}~{}(m=1,\ldots,M).$$

$$\left|{\mathbb{N}}_D^M\right|=\left(\genfrac{}{}{0ex}{}{D+M-1}{D}\right)=\frac{(D+M-1)!}{D!(M-1)!}=O(M^D).$$

$${\mathbb{N}}_D^J:=\left\{(d_1,\dots ,d_M)\in {\mathbb{N}}_D^M\right\}{d}_m=0(m\notin J).$$

$${\mathbf{b}}^J(\mathbf{t}):=\sum _{\mathbf{d}\in {\mathbb{N}}_D^J}\left(\genfrac{}{}{0ex}{}{D}{\mathbf{d}}\right){\mathbf{t}}^{\mathbf{d}}{\mathbf{p}}_{\mathbf{d}}.$$

$$\mathbf{t}\in {\Delta }^J\subseteq {\Delta }^I,$$

$${\mathbf{b}}^I(\mathbf{t})=\sum _{\mathbf{d}\in {\mathbb{N}}_D^I}\left(\genfrac{}{}{0ex}{}{D}{\mathbf{d}}\right){\mathbf{t}}^{\mathbf{d}}{\mathbf{p}}_{\mathbf{d}}=\sum _{\mathbf{d}\in {\mathbb{N}}_D^J}\left(\genfrac{}{}{0ex}{}{D}{\mathbf{d}}\right){\mathbf{t}}^{\mathbf{d}}{\mathbf{p}}_{\mathbf{d}}={\mathbf{b}}^J(\mathbf{t}).$$

$$X_J:=\left\{\mathbf{x}\in X\right\}{\mathbf{f}}_J(y)\nprec {\mathbf{f}}_J(x)\text{for all}\mathbf{y}\in X,$$

Full text

Introduction

A multi-objective optimization problem is a problem that minimizes multiple objective functions $f_{1},\dots,f_{M}:X\to\mathbb{R}$ over a common domain $X\subseteq\mathbb{R}^{L}$:

[TABLE]

Different functions usually have different minimizers, and one needs to consider a trade-off that two solutions $\bm{x},\bm{y}\in X$ may satisfy $f_{i}(\bm{x})<f_{i}(\bm{y})$ and $f_{j}(\bm{x})>f_{j}(\bm{y})$. According to Pareto ordering, i.e.,

[TABLE]

the goal of multi-objective optimization is to obtain the Pareto set

[TABLE]

and the Pareto front

[TABLE]

which describe the best-compromising solutions and their values of the conflicting objective functions, respectively.

In industrial applications, obtaining the whole Pareto set/front rather than a single solution enables us to compare promising alternatives and to explore new innovative designs, whose concept is variously refered to as innovization (?), multi-objective design exploration (?) and design informatics (?). Quite a few real-world problems involve simulations and/or experiments to evaluate solutions (?) and lack the mathematical expression of their objective functions and derivatives. Multi-objective evolutionary algorithms are a tool to solve such problems where the Pareto set/front is approximated by a population, i.e., a finite set of sample points (?).

While the available sample size is limited due to expensive simulations and experiments, it is well-known that the dimensionality of the Pareto set/front increases as the number of objectives grows. Describing high-dimensional Pareto sets/fronts with small samples is one of the key challenges in many-objective optimization today (?).

A considerable number of real-world applications share an interesting structure: their Pareto sets and/or Pareto fronts are often homeomorphic to an $(M-1)$-dimensional simplex. See for example (?; ?; ?; ?). This observation has been theoretically backed up in some cases (?; ?; ?). A recent study (?) pointed out that for all $m\leq M$, each $(m-1)$-dimensional face of such a simplex is the Pareto set of a subproblem optimizing $m$ objective functions of the original problem.

By exploiting this simplex structure, we considers the problem of fitting a hyper-surface to the Pareto front. In statistics, machine learning and related fields, regression problems have been considered in Euclidean space without boundary or at most with coordinate-wise upper/lower boundaries (?; ?). Hyper-surface models developed in such spaces are not suitable for Pareto fronts since a simplex has a non-axis-parallel boundary whose skeleton structure is described by its faces.

This paper proposes a new model and its fitting algorithm for approximating Pareto fronts with the simplex structure. Our contibution can be summarized as follows:

We define a new class of multi-objective optimization problems called the simplicial problem in which the Pareto set/front have the simplex structure discussed above. We propose a Bézier simplex model, which is a generalization of a Bézier curve (?). 2. 2.

We prove that Bézier simplices can approximate the Pareto set/front (as well as the objective map between them) of any simplicial problem with arbitrary accuracy. 3. 3.

We propose a Bézier simplex fitting algorithm. Exploiting the simplex structure of the Pareto set/front of a simplicial problem, this algorithm decomposes a Bézier simplex into low-dimensional simplices and fits each of them, inductively. This approach allows us to reduce the number of parameters to be estimated at a time. 4. 4.

We evaluate the approximation accuracy of the proposed method with synthetic and real-world optimization problems; compared to a conventional response surface model for Pareto fronts, our method exhibits a better boundary approximation while keeping the almost same quality of interior approximation. As a result, a five-objective Pareto front is described with only tens of sample points.

Preliminaries

Let us introduce notations for defining simplicial problems and review an existing method of Bézier curve fitting.

Simplicial Problem

A multi-objective optimization problem is denoted by its objective map $\bm{f}=(f_{1},\dots,f_{M}):X\to\mathbb{R}^{M}$. Let $I:=\set{1,\ldots,M}$ be the index set of objective functions and

[TABLE]

be the standard simplex in $\mathbb{R}^{M}$. For each non-empty subset $J\subseteq I$, we call

[TABLE]

the $J$-face of $\Delta^{M-1}$ and

[TABLE]

the $J$-subproblem of $\bm{f}$. For each $0\leq m\leq M-1$, we call

[TABLE]

the $m$-skeleton of $\Delta^{M-1}$.

The problem class we are interested in is as follows:

Definition 1** **(Simplicial problem)

A problem $\bm{f}:X\to\mathbb{R}^{M}$ is simplicial if there exists a map $\bm{\phi}:\Delta^{M-1}\to X$ such that for each non-empty subset $J\subseteq I$, its restriction $\bm{\phi}|_{\Delta^{J}}:\Delta^{J}\to X$ gives homeomorphisms

[TABLE]

We call such $\bm{\phi}$ and $\bm{f}\circ\bm{\phi}$ a triangulation of the Pareto set $X^{*}(\bm{f})$ and the Pareto front $\bm{f}X^{*}(\bm{f})$, respectively. For each non-empty subset $J\subseteq I$, we call $X^{*}(\bm{f}_{J})$ the $J$-face of $X^{*}(\bm{f})$ and $\bm{f}X^{*}(\bm{f}_{J})$ the $J$-face of $\bm{f}X^{*}(\bm{f})$. For each $0\leq m\leq M-1$, we call

[TABLE]

the $m$-skeleton of $X^{*}(\bm{f})$ and $\bm{f}X^{*}(\bm{f})$, respectively.

By definition, any subproblem of a simplicial problem is again simplicial. As shown in 1b, the Pareto sets forms a simplex. The second condition asserts that $\bm{f}|_{X^{*}(\bm{f})}:X^{*}(\bm{f})\to\mathbb{R}^{M}$ is a $C^{0}$-embedding. This means that the Pareto front of each subproblem is homeomorphic to its Pareto set as shown in 1c. Therefore, the Pareto set/front of an $M$-objective simplicial problem can be identified with a curved $(M-1)$-simplex. We can find its $J$-face by solving the $J$-subproblem.

The above structure appears in a broad range of applications. In operations research, the Pareto set of the facility location problem under the $L^{2}$-norm is shown to be the convex hull of single-objective optima (?). When the optima are in general position, the Pareto set becomes a simplex. Similar observations are also reported under other norms (?). In economics, the Pareto set of the pure exchange economy with $M$ players is known to be homeomorphic to an $(M-1)$-dimensional simplex (?). In hydrology, the two-objective Pareto set of a hydrologic cycle model calibration is observed to be a curve. Its end points are single-objective optima, and the end points correspond to end points of the Pareto front curve (?). In addition, a recent study pointed out that the Pareto set of an $M$-objective convex optimization problem is diffeomorphic to an $(M-1)$-dimensional simplex and that its $(m-1)$-dimensional faces are the Pareto sets of $m$-objective subproblems for all $m\leq M$ (?).

Bézier Curve Fitting

Since the Pareto front of any two-objective simplicial problem is a curve with two end points in $\mathbb{R}^{3}$, the Bézier curve would be a suitable model for describing it.

In $\mathbb{R}^{M}$, the Bézier curve of degree $D$ is a parametric curve, i.e., a map $\bm{b}:[0,1]\to\mathbb{R}^{M}$ determined by $D+1$ control points $\bm{p}_{0},\ldots,\bm{p}_{D}\in\mathbb{R}^{M}$ (?):

[TABLE]

where $\binom{D}{d}$ represents the binomial coefficient. The parameter $t$ moves from $t=0$ to $t=1$, giving a curve $\bm{b}(t)$ with two end points $\bm{b}(0)=\bm{p}_{0}$ and $\bm{b}(1)=\bm{p}_{D}$.

Given sample points $\bm{x}_{1},\ldots,\bm{x}_{N}\in\mathbb{R}^{M}$, a Bézier curve can be fitted by solving the following problem (?):

[TABLE]

where $t_{n}~{}(n=1,\ldots,N)$ and $\bm{p}_{d}~{}(d=0,\ldots,D)$ are variables to be optimized. Notice that $t_{n}~{}(n=1,\ldots,N)$ are introduced to calculate residuals for each sample point. The error function 2 to be minimized represents the sum of squared residuals for sample points.

If one fix all control points $\bm{p}_{d}~{}(d=0,\ldots,D)$, the Bézier curve $\bm{b}(t)$ is determined. The error function 2 can be now separately optimized by minimizing $\left\|\bm{b}(t_{n})-\bm{x}_{n}\right\|^{2}$ with respect to each $t_{n}$. The solution $t_{n}$ is the foot of a perpendicular line from a sample point $\bm{x}_{n}$ to the Bézier curve $\bm{b}([0,1])$, which satisfies

[TABLE]

Since 3 is a nonlinear equation, Newton’s method is used to find the solution $t_{n}$.

If one fix all parameters $t_{n}~{}(n=1,\ldots,N)$, the Bézier curve $\bm{b}(t_{n})$ becomes a linear function with respect to the control points $\bm{p}_{0},\ldots,\bm{p}_{D}$. The error function 2 can be now optimized by solving linear equations with respect to all $\bm{p}_{d}$.

Algorithm 1 shows Borges and Pastva’s method, which alternately adjusts parameters and control points. This algorithm is also used to adjust some of the control points with remaining points fixed (?).

Bézier Simplex Fitting

To describe the Pareto front of an arbitrary-objective simplicial problem, however, Bézier curve fitting is not enough, and we need to generalize it to the Bézier simplex. We propose a method of fitting a Bézier simplex to the Pareto front of a simplicial problem.

Bézier Simplex

Let $\mathbb{N}$ be the set of nonnegative integers and

[TABLE]

For $\bm{t}:=(t_{1},\ldots,t_{M})\in\mathbb{R}^{M}$ and $\bm{d}:=(d_{1},\dots,d_{M})\in\mathbb{N}^{M}$, we denote by $\bm{t}^{\bm{d}}$ a monomial $t^{d_{1}}_{1}t^{d_{2}}_{2}\cdots t^{d_{M}}_{M}$. As shown in Figure 2, the Bézier simplex of degree $D$ in $\mathbb{R}^{M}$ is a map $\Delta^{M-1}\to\mathbb{R}^{M}$ determined by control points $\bm{p}_{\bm{d}}\in\mathbb{R}^{M}$ $(\bm{d}\in\mathbb{N}_{D}^{M})$:

[TABLE]

where $\binom{D}{\bm{d}}$ represents a polynomial coefficient

[TABLE]

Approximation Theorem

Is the Bézier simplex a suitable model for describing the Pareto front of a simplicial problem? Let us check that for any simplicial problem, Bézier simplices can approximate the Pareto front (as well as the Pareto set and the objective map between them) with arbitrary accuracy.

We begin with a more general proposition that any continuous map on a simplex can be approximated by some Bézier simplex as a map:

Theorem 1

Let $\bm{\phi}:\Delta^{M-1}\to\mathbb{R}^{M}$ be a continuous map. There exists an infinite sequence of Bézier simplices $\bm{b}^{(i)}:\Delta^{M-1}\to\mathbb{R}^{M}$ such that

[TABLE]

The proof of Theorem 1 is shown in Appendix A. Recall Definition 1 that a simplicial problem $\bm{f}:X\to\mathbb{R}^{M}$ admits a triangulation of the Pareto set $\bm{\phi}:\Delta^{M-1}\to X^{*}(\bm{f})$, which induces a triangulation of the Pareto front $\bm{f}\circ\bm{\phi}:\Delta^{M-1}\to\bm{f}X^{*}(\bm{f})$. These triangulations can be $\bm{\phi}$ in Theorem 1. Addtionally, the restricted objective map $\bm{f}:X^{*}(\bm{f})\to\bm{f}X^{*}(\bm{f})$ has the graph

[TABLE]

and $\bm{\phi}:\Delta^{M-1}\to X^{*}(\bm{f})$ also induces a triangulation of the graph

[TABLE]

since $\bm{f}:X^{*}(\bm{f})\to\bm{f}X^{*}(\bm{f})$ is continuous by definition and any continuous map induces a homeomorphism from the domain to the graph. We thus get the desired result.

Corollary 2

Let $X^{*}$ be the Pareto set, the Pareto front or the graph of the objective map restricted to the Pareto set of a simplicial problem. There exists an infinite sequence of Bézier simplices $\bm{b}^{(i)}:\Delta^{M-1}\to\mathbb{R}^{M}$ such that

[TABLE]

where $d_{\mathrm{H}}$ is the Hausdorff distance and $B^{(i)}$ are images of Bézier simplices: $B^{(i)}:=\bm{b}^{(i)}(\Delta^{M-1})$.

All-at-Once Fitting

Let us consider algorithms to realize such a sequence of Bézier simplices. First, we develop a straightforward generalization of the Bézier curve fitting method. We call this method the all-at-once fitting.

Given sample points $\bm{x}_{1},\ldots,\bm{x}_{N}\in\mathbb{R}^{M}$, a Bézier simplex can be fitted by solving the following problem, which is a multi-dimensional analogue of the problem 2:

[TABLE]

where $\bm{t}_{n}=(t_{n1},\ldots,t_{nM})~{}(n=1,\ldots,N)$ and $\bm{p}_{\bm{d}}~{}(\bm{d}\in\mathbb{N}^{M}_{D})$ are variables to be optimized.

As is the case of Bézier curve fitting, if one fix all control points $\bm{p}_{\bm{d}}~{}(\bm{d}\in\mathbb{N}^{M}_{D})$, the Bézier simplex $\bm{b}(\bm{t})$ is determined. The error function 5 can be now separately optimized by minimizing $\left\|\bm{b}(\bm{t}_{n})-\bm{x}_{n}\right\|^{2}$ with respect to each $\bm{t}_{n}$. The solution $\bm{t}_{n}$ is the foot of a perpendicular line from a sample point $\bm{x}_{n}$ to the Bézier simplex $\bm{b}(\Delta^{M-1})$, which satisfies

[TABLE]

Since 6 is a system of nonlinear equations, Newton’s method is used to find the solution $\bm{t}_{n}$.

If one fix all parameters $\bm{t}_{n}~{}(n=1,\ldots,N)$, the Bézier simplex $\bm{b}(\bm{t}_{n})$ is a linear function with respect to the control points $\bm{p}_{\bm{d}}~{}(\bm{d}\in\mathbb{N}^{M}_{D})$. The error function 5 can be now optimized by solving linear equations with respect to all $\bm{p}_{\bm{d}}~{}(\bm{d}\in\mathbb{N}^{M}_{D})$.

As well as Borges and Pastva’s method, the all-at-once fitting alternatively adjusts parameters and control points. We describe the all-at-once fitting in Algorithm 2. This algorithm is also used to adjust some of the control points with remaining control points fixed.

While a Bézier simplex has enough flexibility as described in Theorem 1, the required sample size for the all-at-once fitting grows quickly with respect to $M$. Let us consider the case of fitting a Bézier simplex of degree $D$ to the Pareto front of an $M$-objective problem. In this case, the number of control points to be estimated is

[TABLE]

Practically, the degree of a Bézier simplex to be fitted is often set as $D=3$; nevertheless the number of control points to be estimated at a time becomes unreasonable for large $M$.

Inductive Skeleton Fitting

To reduce the required sample size, we consider decomposing a Bézier simplex into subsimlices and fitting each subsimplex one by one from low dimension to high dimension. This approach allows us to reduce the number of control points to be estimated at a time. We call this method the inductive skeleton fitting.

Let $I:=\Set{1,\ldots,M}$ and for each non-empty subset $J\subseteq I$, we define

[TABLE]

The $J$-face of an $(M-1)$-Bézier simplex of degree $D$ is a map $\Delta^{J}\to\mathbb{R}^{M}$ determined by control points $\bm{p}_{\bm{d}}\in\mathbb{R}^{M}$ for all $\bm{d}\in\mathbb{N}_{D}^{J}$:

[TABLE]

For arbitrary parameters $\bm{t}$ satisfying

[TABLE]

all entries that are not in $J$ are 0. It then holds that

[TABLE]

This means that, for each non-empty subset $J\subset I$, the $J$-face of the $(M-1)$-Bézier simplex of degree $D$ is the $(\left|J\right|-1)$-Bézier simplex of the same degree and it is determined by the control points $\bm{p}_{\bm{d}}$ of the $(M-1)$-Bézier simplex satisfying $\bm{d}\in\mathbb{N}_{D}^{J}$. To exploiting this structure, the inductive skeleton fitting decomposes a Bézier simplex into subsimplices and fits a subsimplex $\bm{b}^{J}(\bm{t})$ to $X_{J}$ for each non-empty subset $J\subseteq I$ in ascending order of cardinality of $J$.

When a problem $\bm{f}=(f_{1},\ldots,f_{M}):X\to\mathbb{R}^{M}$ is simplicial, we can provide a set of subsamples for running the inductive skeleton fitting. Remember that the Pareto front $\bm{f}X^{*}(\bm{f})$ has the same skeleton as the $(M-1)$-simplex: $\bm{f}X^{*}(\bm{f}_{J})\subseteq\bm{f}X^{*}(\bm{f}_{I})$ for all $\emptyset\neq J\subseteq I$. Thus given a sample $X$ of $\bm{f}X^{*}(\bm{f})$, we can decompose it into subsamples

[TABLE]

each of which represents a sample of the $J$-face $\bm{f}X^{*}(\bm{f}_{J})~{}(\emptyset\neq J\subseteq I)$. Algorithm 3 summarizes the inductive skeleton fitting.

Unlike the all-at-once fitting, the inductive skeleton fitting allows us to reduce the number of control points to be estimated at a time. Let us consider the case of fitting a subsimplex $\bm{b}^{J}(\bm{t})$ with $\left|J\right|=m$. As we described before, the subsimplex $\bm{b}^{J}(\bm{t})$ has $\left|\mathbb{N}^{m}_{D}\right|=\binom{D+m-1}{D}=\frac{(D+m-1)!}{D!(m-1)!}$ control points. In practice, it is sufficient to set $D$ as a small value compared to $M$. In such a case, the inductive skeleton fitting estimates at most $D$-objective solutions, then we only have to adjust at most $\max_{m=1,\ldots,D}\binom{D+m-1}{D}$ control points for each step. Notice that this number does not depend on $M$ but on $D$. Therefore, in case of the inductive skeleton fitting, the number of control points to be estimated at a time is much smaller than $\binom{D+M-1}{D}$ for large $M$.

Numerical Experiments

The small sample behavior of the proposed method is examined using Pareto front samples of varying size.111The source code is available from https://github.com/rafcc/stratification-learning.

Data Sets

To investigate the effect of the Pareto front shape and dimensionality, we employed six synthetic problems with known Pareto fronts, all of which are simplicial problems. Schaffer, ConstrEx and Osyczka2 are two-objective problems. Their Pareto fronts are a curved line that can be triangulated into two vertices and one edge. 3-MED and Viennet2 are three-objective problems. Their Pareto fronts are a curved triangle that can be triangulated into three vertices, three edges and one face. 5-MED is a five-objective problem. Its Pareto front is a curved pentachoron that can be triangulated into five vertices, ten edges, ten faces, five three-dimensional faces and one four-dimensional face. We generated Pareto front samples of 3-MED and 5-MED by AWA(objective) using default hyper-parameters (?). Pareto front samples of the other problems were taken from jMetal 5.2 (?).

To assess the practicality of the proposed method, we also used one real-world problem called S3TD. This is a four-objective problem222The original problem is five-objective, but we dropped one objective since there are two objectives that are highly correlated. of designing a silent super-sonic aircraft. The Pareto front has not been exactly known, and we only have an inaccurate sample of 58 points obtained in the prior study (?). Its simpliciality is also unknown. If it is assumed to be simplicial, then the Pareto front is a curved tetrahedron that can be triangulated into four vertices, six edges, four faces, one three-dimensional face. The problem definitions of the synthetic and real-world problems are described in Appendix B.

For each problem, the Pareto front sample is split into a training set and a validation set. The training set is further decomposed into subsamples for the Pareto fronts of subproblems: each $m$-objective subproblem has a subsample consisting of $N_{m}$ points. Experiments were conducted on all combinations of the following subsample sizes:333We fix $N_{4}=N_{5}=0$ since all methods used in our experiments do not require these subsamples to fit the Pareto front.

[TABLE]

An $M$-objective problem has $M$ single-objective problems, $M(M-1)/2$ two-objective problems and $M(M-1)(M-2)/6$ three-objective problems. The total sample size of the training set is

[TABLE]

The details of making Pareto front subsamples are described in Appendix C.

Methods

We compared the following surface-fitting methods:

• •

the inductive skeleton fitting (Algorithm 3);

• •

the all-at-once fitting (Algorithm 2);

• •

the response surface method (?).

The inductive skeleton fitting used each training subsample for each subproblem. The all-at-once fitting and the response surface method used all training subsamples as a whole. We set the initial control points $\bm{p}_{\bm{d}}^{0}$ $(\bm{d}\in\mathbb{N}_{D}^{J},\ J\subset I,\ \left|J\right|=1)$, which are the vertices of the Bézier simplex, to be the single-objective optima, and the rest of control points were set to be the simplex grid spanned by them. Newton’s method used in the inductive skeleton fitting and the all-at-once fitting employed the first and second analytical derivatives, and was terminated when the number of iterations reached 100 or when the following condition is satisfied:

[TABLE]

The fitting algorithm was terminated when the number of iterations reached 100 or when the following condition is satisfied:

[TABLE]

where $\mathrm{SSR}_{i}$ is the value of the loss function 5 at the $i$-th iteration.

According to a prior study (?), the response surface method treated multi-dimensional objective values as a function from the first $M-1$ objective values to the last objective value by using multi-polynomial regression with cubic polynomials. We removed cross terms of degree three so that the number of explanatory variables did not exceed the sample size.

Performance Measures

To evaluate how accurately an estimated hyper-surface approximates the Pareto front, we used the generational distance (GD) (?) and the inverted generational distance (IGD) (?):

[TABLE]

where $X$ is a finite set of points sampled from an estimated hyper-surface and $Y$ is a validation set.

Figure 3 depicts what GD and IGD assess. These values can be viewed as the avarage length of arrows in each plot. 3a implies that GD becomes high when the estimated front has a false positive area which is far from the Pareto front. Conversely, 3b tells that IGD becomes high when the Pareto front has a false negative area which is not covered by the estimated front. Thus, we can say that the estimated hyper-surface is close to the Pareto front if and only if both GD and IGD are small.

For the all-at-once fitting and the inductive skeleton fitting, we sampled the estimated Bézier simplex as follows:

[TABLE]

For the response surface method, we sampled the estimated response surface as follows:

[TABLE]

where $r(x_{1},\ldots,x_{M-1})$ is an estimated polynomial function. We repeated experiments 20 times444We set the number of trials on the basis of the power analysis. According to our preliminary experiments, the effect size of GD/IGD values was approximately $d=0.6$. In this case, the sample size $n=20$ is required to assume the power $>0.8$ with significance level $p=0.05$. with different training sets and computed the average and the standard deviations of their GDs and IGDs.

Results

For each problem and method, the average and the standard deviation of the GD and IGD are shown in Table 1. In Table 1, we highlighted the best score of GD and IGD out of all methods for each problem and added the results of Mann-Whitney’s one-tail U-test to check the best score was smaller than that of the other methods significantly. In case of conducting multiple tests, we corrected p-value by Holm’s method.

For Bézier simplex fitting methods, the number of iterations until termination is also shown. The table shows that the fitting was finished in approximately three iterations in all cases except for the all-at-once fitting on S3TD.

Inductive skeleton vs. all-at-once

For all the two-objective problems, both methods obtained almost identical GD and IGD values. For three- and five-objective problems, the inductive skeleton fitting significantly outperformed the all-at-once fitting in both GD and IGD. Especially in Viennet2, the all-at-once fitting exhibited unstable behavior in GD while the inductive skeleton fitting was stable.

Inductive skeleton vs. response surface

The inductive skeleton fitting achieved better IGDs in Schaffer, Osyczka2, 3-MED, and 5-MED but worse values in ConstrEx and Viennet2. In terms of the average GD, the inductive skeleton fitting was better in all the problems except ConstrEx.

Discussion

In this section, we discuss approximation accuracy, required sample size, practicality for real-world data and objective map approximation of our method.

Approximation Accuracy

For three- and five-objective problems, the inductive skeleton fitting achieved slightly better IGDs. The inductive skeleton fitting adjusts a small subset of control points at each step with already-adjusted control points of the faces. This reduces the number of parameters estimated at a time, which seems to prevent over-fitting.

More significant differences can be seen in GD. Figure 4 shows the cause of these defferences: The all-at-once fitting obtained an overly-spreading Bézier simplex while the inductive skeleton fitting found an exactly-spreading Bézier simplex. Minimizing the squared loss 5 only imposes that all sample points are close to the Bézier simplex, which leads to a good IGD. However, it does not impose that all Bézier simplex points are close to the sample points, which is needed to achieve a good GD. The inductive skeleton fitting stipulates that each face of the Bézier simplex must be close to each face of the Pareto front (i.e., the Pareto front of each subproblem), which minimizes GD.

Similarly to the all-at-once fitting, the response surface has poor GDs. As we do not know the true Pareto set, it is difficult for grid sampling to obtain an exactly-spreading surface.

ConstrEx and Viennet2 are the only problems where the inductive skeleton fitting was worse than the response surface method. ConstrEx is a non-smooth curve that cannot be fitted by a single Bézier curve (see Appendix D in Appendix D). This type of Pareto fronts leads to a challenging problem: developing a method for gluing multiple Bézier simplices to express a non-smooth surface. Viennet2 is a smooth surface but its curvature is severely sharp (see Appendix D in Appendix D). Although Viennet2 is defined by quadratic functions (see Appendix B), its Pareto front cannot be fitted by the Bézier simplex of degree three. This means that setting the degree of the Bézier simplex grater than or equal to the maximum degree in the problem definition does not ensure approximation accuracy. It would be useful to develop a way to understand the required degree of the Bézier simplex from the problem definition.

Required Sample Size

To achieve good accuracy, how many sample points does the Bézier simplex require? Figure 6 shows transitions of GD and IGD on 5-MED when the sample size of each three-objective subproblem (two-dimensional face) varies $N_{3}=1,\dots,10$ with fixed $N_{1}=1$ and $N_{2}=2$. Surprisingly, both GD and IGD have already converged at $N_{3}=4$.

In case of $D=3$ for 5-MED, which is a five-objective ($M=5$) problem, the number of all control points to be esitimated is $\binom{3+5-1}{5}=35$. While the all-at-once fitting fits them simultaneously, in case of the induuctive skeleton method, we only have to fit two and one control points for each two- and three-objective subproblem respectively.This result demonstrates that the inductive skeleton fitting allows us to fit Pareto fronts with small samples by reducing the number of control points to be estimated at a time.

Practicality for Real-World Data

Now we discuss the practicality of the inductive skeleton fitting. Firstly, we focus on the performance of S3TD, which is a real-world problem. Table 1 indicates that the average IGD of the inductive skeleton fitting on S3TD was worse than that of the response surface method. However, the IGD itself is an order of $10^{-1}$, which small in absolute sense considering the number of the objective functions ($M=4$) and the sample size of the validation set ($N=58$) on S3TD. The average GD of the inductive skeleton method was smaller than that of the other methods and the difference was significant with significance level $p=0.1$. Figure 6 shows transitions of GD and IGD on S3TD when the sample size of each three-objective subproblem (two-dimensional face) varies $N_{3}=1,\dots,10$ with fixed $N_{1}=1$ and $N_{2}=5$. We can observe that both GD and IGD already converged at $N_{3}=4$ as Figure 6. These results suggest that the inductive skeleton fitting can work properly in real-world problems.

In real-world problems, it seems to be difficult to know that the Pareto set/front is a simplex and in some cases, it is degenerated. However, a statistical test has been proposed to judge whether the Pareto set/front of a given problem is a simplex or not (?). This test tells us whether we can apply the inductive skeleton fitting to real-world problems. Furthermore, even if the Pareto front is degenerated, Theorem 1 states the Bezier simplex model can fit degenerate Pareto sets/fronts.

Objective Map Approximation

As well as the Pareto front approximation discussed and demonstrated so far, Theorem 1 ensures that our method can be used for approximating the Pareto set, $X^{*}(\bm{f})$, and the restricted objective map, $\bm{f}:X^{*}(\bm{f})\to\bm{f}X^{*}(\bm{f})$. Such approximations may provide richer information about the problem. To check the performance of graph approximation, we made a sample of solution-objective pairs of 5-MED $X=\Set{(\bm{x},\bm{f}(\bm{x}))\in\mathbb{R}^{5}\times\mathbb{R}^{5}}{\bm{x}\in X^{*}(\bm{f})}$ and applied the inductive skeleton fitting to it. Other settings were the same as the Pareto front approximation. Table 2 shows the GD and IGD values where the results of the all-at-once fitting are just for baseline. As the inductive skeleton fitting achieved GD and IGD of $10^{-1}$ order of magnitude, which means that it accurately approximates the graph of the restricted objective map in the abslute sense. In 5-MED, the graph of the restricted objective map on the Pareto set is a four-dimensional topological simplex in a 10-dimensional space. Although the codimension is six times higher than the case of Pareto front approximation, we have similar GD and IGD to ones of 5-MED in Table 1. This fact implies that the approximation error mainly depends on the intrinsic dimensionality rather than the ambient dimensionality.

Conclusions

In this paper, we have proposed the Bézier simplex model and its fitting algorithm for approximating Pareto fronts of multi-objective optimization problems. An approximation theorem of the Bézier simplex has been proven. Numerical experiments have shown that our fitting algorithm, the inductive skeleton fitting, obtains exactly-spreading Bézier simplices over synthetic Pareto fronts of different shape and dimensionality. It has been also observed that a real-world problem with four objectives but only 58 points is accurately fitted. The proposed model and its fitting algorithm drastically reduce the sample size required to describe the entire Pareto front.

The current algorithm has some drawbacks. The inductive skeleton fitting requires that each subproblem must have a non-empty sample, which may be too demanding for some practice. Furthermore, the loss function for fitting does not take into account of the sampling error of the Pareto front. To get a more stable approximation under wild samples, we plan to extend the (deterministic) Bezier simplex to a probabilistic model.

We believe that the potential applications of the method are not limited to the use in post-optimization process. It will be introduced in evolutionary algorithms and dynamically evolved to accelerate search. We also expect that this model can be applied to regression problems that have complex boundary conditions such as learning from multi-labeled data.

Acknowledgements

We wish to thank Dr. Yuichi Ike for checking the proof and making a number of valuable suggestions.

Appendix A Appendix A: Proof of Theorem 1

In this section, we prove that Bézier simplices can approximate pareto sets of simplicial problems in any accuracy. Let $\bm{b}:\Delta^{M-1}\to\mathbb{R}^{M}$ be a Bézier simplex.We say this function $\bm{b}:\Delta^{M-1}\to\mathbb{R}^{M}$ Bézier maps and the projections to its coordinates $pr\circ\bm{b}:\Delta^{M-1}\to\mathbb{R}$ Bézier functions. Namely,

[TABLE]

where $b_{1},\ldots,b_{M}$ are Bézier functions.

**Theorem **1

Let $\bm{\phi}:\Delta^{M-1}\to\mathbb{R}^{M}$ be a $C^{0}$-embedding. There exists an infinite sequence of Bézier simplices $\bm{b}^{(i)}:\Delta^{M-1}\to\mathbb{R}^{M}$ such that

[TABLE]

Remark 1

By definition, if $X^{*}$ is the Pareto set of a simplicial problem, then there is a homeomorphism $\bm{\phi}:\Delta^{M-1}\to X^{*}$. Therefore, Theorem 1 implies that Pareto sets of simplicial problems can be approximated by Bézier simplices.

To prove this theorem, we use the following theorem by Stone-Weierstrass. Let us prepare some notations before stating Stone-Weierstrass theorem. Let $X$ be a compact metric space. Let $C(X)$ be the set of continuous functions over $X$. We regard $C(X)$ as a Banach algebra with the sup-norm,

[TABLE]

A subalgebra of $C(X)$ is a vector subspace $S$ such that for any $f,g\in S$, $f\cdot g\in S$. We say that a subalgebra $S$ separates points if for any different points $\bm{x},\bm{y}\in X,$ there exists $f\in S$ so that $f(\bm{x})\neq f(\bm{y})$. We say that $S$ is unital if $S$ contains $1$.

Theorem 3** **(Stone-Weierstrass)

Let $S$ be a unital subalgebra of $C(X)$. Assume $S$ separates any two points in $X$. Then $S$ is dense in $C(X)$, namely $\bar{S}=C(X)$.

We use Stone-Weierstrass theorem in the case where $X=\Delta^{M-1}$ and $S$ is the set of Bézier functions. To confirm the assumptions of Stone-Weierstrass theorem, we need the following propositions.

Proposition 1

Let $S$ be the set of Bézier functions of any degree. Then $S$ is a unital subalgebra of $C(\Delta^{M-1})$.

**Proof **

We can easily see that the product of Bézier functions of degree $D_{1}$ and $D_{2}$ is a Bézier function of degree $D_{1}+D_{2}$. Hence, it is enough to show that the linear combination of Bézier functions is a Bézier function again. Since it is clear that the linear combination of Bézier functions of the same degree is a Bézier function, it is enough to show that Bézier functions of degree $D_{1}$ can be represented by Bézier functions of degree $D_{2}$, if $D_{1}\leq D_{2}$. To see this, let us denote by $S_{D}$ the set of the Bézier functions whose degree is equal to $D$. We can regard $S_{D}$ as a linear subspace of $\mathbb{R}[t_{1},\ldots,t_{M}]$. Let $S^{\prime}_{D}$ be the set of polynomials of degree lower than or equal to $D$ in $\mathbb{R}[t_{1},\ldots,t_{M}]$ (the degree is the total degree by $t_{1},\ldots,t_{M}$). By the definition of Bézier functions, $S_{D}$ is contained in $S^{\prime}_{D}$. Since $t_{1}^{D_{1}},\ldots,t_{M}^{D_{M}}$ are linearly independent over $\mathbb{R}$, the dimension of $S_{D}$ is equal to $\binom{M+D-1}{M-1}$. However, the dimension of $S^{\prime}_{D}$ is equal to $\binom{M+D-1}{M-1}$. Therefore, $S_{D}$ is equal to $S^{\prime}_{D}$. Hence, $S_{D_{2}}=S^{\prime}_{D_{2}}\supset S^{\prime}_{D_{1}}=S_{D_{1}}$. $S$ is a subalgebra of $C(\Delta^{M-1})$. Clearly $S^{\prime}_{0}$ contains $1$. This implies that $S$ is unital.

Proposition 2

Let $S$ be the set of Bézier functions of any degree. Then $S$ separates points.

**Proof **

Take two different points $\bm{s},\bm{t}\in\Delta^{M-1}$. Consider the Bézier functions of degree $1$. By the definition of Bézier functions, Bézier functions of degree $1$ define arbitrary hyperplanes (or equivalently linear equations). We can find a hyperplane $\bm{b}$ which contains $\bm{s}$ and does not contain $\bm{t}$. In this case, $\bm{b}(\bm{s})\neq 0=\bm{b}(\bm{t})$. Hence $S$ separates points.

Now we prove the main theorem.

**Proof **** **(proof of Theorem 1)

The Hausdorff distance is defined as

[TABLE]

Hence, if there is a Bézier simplex $\bm{b}:\Delta^{M-1}\to\mathbb{R}^{M}$ such that for any $\bm{t}\in\Delta^{M-1}$, $\left|\bm{\phi}(\bm{t})-\bm{b}(\bm{t})\right|$ is enough small, we obtain the assertion. Let $S$ be the set of Bézier functions of arbitrary degree. By Propositions 1 and 2, $S$ is a unitary algebra which separates points on $\Delta^{M-1}$. By Theorem 3, $S$ is dense in $C(\Delta^{M-1})$. Hence we can find Bézier functions $b_{i}$ such that $\left\|pr_{i}\circ\bm{\phi}-b_{i}\right\|_{\infty}$ is enough small, where $pr_{i}$ denote the projection to the $i$-th coordinate. Put $\bm{b}:=(b_{1},\ldots,b_{M})$. Then

[TABLE]

holds. If the right hand side is enough small, $\sup_{\bm{t}\in\Delta^{M-1}}\left|\bm{\phi}(\bm{t})-\bm{b}(\bm{t})\right|$ is also enough small.

Appendix B Appendix B: Problem Definition

Schaffer

is a one-variable two-objective problem defined by:

[TABLE]

ConstrEx

is a two-variable two-objective problem defined by:

[TABLE]

Osyczka2

is a six-variable two-objective problem defined by:

[TABLE]

Viennet2

is a two-variable three-objective problem defined by:

[TABLE]

$M$-MED

is an $M$-variable $M$-objective problem defined by:

[TABLE]

S3TD

is a 58-variable four-objective problem of designing a super-sonic airplane. The original problem has five objectives:

Minimizing pressure drug 2. 2.

Minimizing fliction drug 3. 3.

Maximizing lift 4. 4.

Minimizing sonic boom intensity 5. 5.

Minimizing structural weight

All objectives are converted into minimization. We removed $f_{1}$ since $f_{1}$ and $f_{4}$ are strongly correlated and $f_{4}$ is the primary objective in this project.

Appendix C Appendix C: Data Description

The description of data that we used in the numerical experiments is summarized in Table 3: Table 3 shows the sample size of each subproblem. To make a data set for each problem, we chose mutually exclusive subsamples of the Pareto fronts of the subproblems.

Appendix D Appendix D: All Results

The Estimated Bézier Simplices and Response Surfaces for Pareto Front

Appendix D shows pairplots of the estimated Bézier simplices and response surfaces for two-objective problems with sample size $(1,3)$. For three- and five-objective problems, the results of 3-MED, Viennet2 and 5-MED with sample size $(1,2,2)$ are shown in Appendices D, D and D, respectively. For a real-world problem, the results of S3TD, which is four-objective problem with sample size $(1,5,5)$ is shown in Appendix D.

Transitions of GD and IGD

Appendix D shows boxplots of GD and IGD when we changed the sample size on the edges for two-objective problems. The results of 3-MED, Viennet2 and 5-MED are shown in Appendices D, D and D, respectively.

Figure 7: Bézier curves and a response surface for Schaffer (left), ConstrEx (center) and Osyczka2 (right) with sample size $(1,3)$.

Figure 8: Bézier triangles and a response surface for 3-MED with sample size $(1,2,2)$.

Figure 9: Bézier triangles and a response surface for Viennet2 with sample size $(1,2,2)$.

Figure 10: Bézier triangles and a response surface for 5-MED with sample size $(1,2,2)$.

Figure 11: Bézier triangles and a response surface for S3TD with sample size $(1,5,5)$.

Figure 12: Sample size $N_{2}$ vs. GD/IGD (boxplot over ten trials).

Figure 13: Sample size $N_{3}$ vs. GD/IGD on 3-MED (boxplot over ten trials).

Figure 14: Sample size $N_{3}$ vs. GD/IGD on Viennet2 (boxplot over ten trials).

Figure 15: Sample size $N_{3}$ vs. GD/IGD on 5-MED (boxplot over ten trials).