Bézier Flow: a Surface-wise Gradient Descent Method for Multi-objective Optimization
Authors: Akiyoshi Sannai
Yasunari Hikima
Ken Kobayashi
Akinori Tanaka
Naoki Hamada
Venue: Accepted by TMLR ()

Abstract: This paper proposes a framework to construct a multi-objective optimization algorithm from a single-objective optimization algorithm by using the Bézier simplex model. Additionally, we extend the stability of optimization algorithms in the sense of Probably Approximately Correct (PAC) learning and define the PAC stability. We prove that it leads to an upper bound on the generalization error with high probability.
Furthermore, we show that multi-objective optimization algorithms derived from a gradient descent-based single-objective optimization algorithm are PAC stable. We conducted numerical experiments with synthetic and real multi-objective optimization problem instances and demonstrated that our method achieved lower generalization errors than the existing multi-objective optimization algorithms.

Published in Transactions on Machine Learning Research (04/2025)
Bézier Flow: a Surface-wise Gradient Descent Method for Multi-objective Optimization
Akioosii Sannai
sannai.ckigoshi. htw.:idddiipshh.ss@htoto....o
Kyoto University
National Institute of Informatics
Yasunari Hikima
hukima. hikkima.yysuunarri@uuiuttuuu com
Fujilsu Limited
Kyusha University
Ken Kobayashi Institute of Science Tokyo
kobayashi.k eng. isc.issct..o ac.jp
Akinori Tanaka
Aksnoori.aaaate @riken..j
RIKEN AIP
RIKEN 3THEMS
Keio University
Naoki Hamada
hamada-n&kkab. com
KLab Inc.
Reviewed on OpenReview: https: openreview. net/ forum? id- IigALubRaj
Abstract
This paper proposes a framework to construct multi-objective optimization algorithm from singlo-objective optimization algorithm by using the Bézier simplex model. Additionally, we exlend the stability of optimization algorithms in the sense of Probably Approximately Correct (PAC) learning and define the PAC stability. We prove that it leads to an upper bound on the generalization error with high probability. Furthermore, we show that multi- objective optimization algorithms derived from gradient descent-based single-objective optimization algorithm PAC stable. We conducted numerical experiments with synthetic and real multi-objective optimization problem instances and demonstrated that method achieved lower generalization errors than the existing multi-objective optimization algorithms.
1 Introduction
A multi-objective optimization problem a problem to seck a solution which minimizes (or maximizes) multiple objective functions f1...., fM +R
R simultancously over a domain X CRR
minimiie
= (())=)))))...,(
subject to
x E X
Each objective function can have different optimal solution, so we need to consider the trade-off between two or more solutions. Therefore, the notion of Parcto ordering is taken into consideration, which is defined by
f(x) f(y) def fm( S (v) for all m 1..... M, and fx) fm (3) for some 1, ... M.

---

Published in Transactions on Machine Learning Research (04/2025)
24
18
(3)
10 (x
(1.a)
(2,3)
(((aa))))
homeomorphism
fo xx ((y))
(1,3.))
(2)
((A))==
X X"
(f
((())
embedding
(f
5(x (
12
Aq1,2)
fi
r(x (1))))))))
(a) Simplex A?.
(b) Pareto set (f).
(c) Pareto front (f)..
Figure 1: A simplicial problem f= (f1, fa) R3 > An M-objective problem I simplicial if
the following conditions are satisficd: (i) there exists homcomorphism $ AM-1 X*) such that
((r) X (ffr for all I C (1,..., M]; (ii) the restriction flx-(n) (J) - RM is topological embedding
(and thus is f o. A AM-1_ - I)..
In multi-objective optimization, the goal is to obtain the Pareto set and Pareto front, which are respectively defined as:
x²f)= E fty A f(x) for all y E X). fX*() (f())
xxEx)()))
The Pareto set set/front usually has an infinite number of points, whereas most of the numerical methods for solving the problem give us a finite set of points as an approximation of the Parcto set set/ set/front (c.g., goal programming (Eichfelder, 2008; Miettinen, 1999), evolutionary computation (Deb, 2001; Deb & Jain, 2013; Zhang& & Li, 2007), homotopy methods (Harada et al., 2007; Hillermeier, 2001), and Bayesian optimization (Hernanezz-Loatoo et al., 2016; Yang et al., 2019)). Such a finite-point approximation cannot reveal the complete shape of the Pareto set and Pareto front. In addition, the finite-point approximation suffers from the "curse of dimensionality" since the dimensionality of the Parelo set and Pareto front. is M - in generic problems (see Wan (1977; 1978) for a rigorous statement). With this background, we focus on an optimization algorithm to obtain a parametric hypersurface describing the Parcto set.
There is a common structure of the Pareto / across wide variety of problems, which can be utilized to enhance approximation. In many problems, obtained solutions imply the Pareto set /front is curved (M - 1) -)simp,,,,, c.g., airplane design (Mastroddi & Gemma, 2013), hydrologic modeling (Vrugt et al., 2003), PI controller tuning (Reynoso-Meza et al., 2015), building design (Safarzadegan Gilan et al., 2016), motor design (Contreras et al., 2016), and Lasso's syper-parameter tuning (Hamada & Ichiki, 2020). To mathematically identify such a class of problems, Kobayashi et al. (2019) defined the simplicial problem (see Figure 1). Hamada et al. (2020) showed that strongly convex problems are simplicial under mild conditions, which implies that facility location (Kuhn, 1967) and phenotypic divergence modeling in evolutionary biology (Shoval et al., 2012) are simplicial. Kobayashi et al. (2019) showed that the Pareto set and Pareto front of any simplicial problem can be approximated with arbitrary accuracy by a Bézier simplex.
By using this advantage of the Bézier simplex model, we propose novel strategy to construct a multi- objective optimization algorithm from a single-objective optimization algorithm. With a given single objective oplimization algorithm, such as a gradient descent method, this scheme updates the Bézier simplex to obtain the Pareto set. In addition, we analyze the theoretical properties of the multi-objective optimization algorithm derived from our scheme. Specifically, we define Probably Approximately Correct (PAC) stability as an extension of the stability of optimization algorithms and prove that the PAC stability leads to an upper bound on the generalization gap in the sense of PAC learning. Our contributions are summarized as follows:
1. We devise a strategy to construct a multi-objective optimization algorithm from a single-objective optimization algorithm with the Bézier simplex. Unlike most of the existing multi-objective optimiza- tion methods, the algorithm derived from our scheme has the advantage of obtaining a parametric hypersurface that represents the Pareto set of a given simplicial, Lipschitz continuous, differentiable multi-objective optimization problem to be solved.
2

---

Published in Transactions on Machine Learning Research (04/2025)
2. We define PAC stability, which is an extension of the stability introduced by Hardt et al. (2016) to the PAC learning settings, and show that PAC stability gives an upper bound on the generalization gap with a high probability. Also, we prove that when we employ a gradient-based optimization algorithm as a single optimization algorithm, the derived multi-objective optimization algorithm is PAC stable.
3. We conducted numerical experiments and demonstrated that the multi-objective optimization algorithm constructed by our scheme achieved lower generalization errors than the existing multi- objective optimization algorithm. In addition, the algorithm given by our scheme can efficiontly obtain the Parelo set with small number of sample points.
Related Work Kobayashi et al. (2019) proposed Bèzier simplex litting algorithms, the all-at-once fitting, and inductive skeleton fitting to describe Pareto fronts, and Tanaka et al. (2020) analyzed the asymptotic risk of the fitting algorithms. The two fitting algorithms focus on post-optimization processes and assume that we have an approximate solution set of the Pareto set in advance. Thus, these algorithms alone cannot solve multi-objective optimization problems. Recently, Maree el. al. (2020) proposed bi-objective optimization algorithm that updates the Bézier curve. However, this algorithm exploits the structure of the bi-objective optimization problem and can not be applied when the number of objective functions is more than or equal to three. To the best of our knowledge, we are the first to propose general framework of multi-objectivee optimization with the Bézier simplex and show its theoretical properties.
2 Preliminaries
2.1 Probability simplex
Let [M] = f1.....,M M) be set of M points. We consider the set of probability distribution t E M over [M]. The set of probability distributions over [M] equal to the simplex
A M-1 AM-1=
-
{ (
(ot.....t
(M) E RM
tm
Σ tm m
M
1m
Let C(X) be the space of continuous functions over X, and we define the function F: [M] + C(X) by F(m) = fm- Then, we have the expectation function
B(F):
C(X)
EF(F)
Furthermore, if fm is strongly convex for all E [M], then the following function is well-defined:
argmin B(f): A
X
W
argmin (F
Note that EF(F) Em fmfm follows from the definition. (F) corresponds to the sum of a function chosen continuously along from f. As a direct consequence from Theorem 2 in Mizota et al. (2021), the mapping argmin E(f) gives a continuous surjection onto X*(f) if fm is strongly convex for all m [M].
2.2 Simplicial problem
A multi-objective optimization problem is characterized by its objective map f = .fi,..., f.)..f))TT x X
We define the J J-subsimpler for an index set S [M] by A ((t....., (T) E A M-1 Imm =0 (m E J)). The problem class we wish consider is problem in which the Pareto set set/front /front has the simplex structure. Such problem class defined as follows.
3

---

Published in Transactions on Machine Learning Research (04/2025)
Definition 2.1 (Kobayashi et al. (2019)). A problem f: X + RM is simplicial if there exists map Φ: AM-1 - X such that for each non-emply subset J = [M], its restriction Φ A ( A - X gives
homeomorphisms
fo AJ + fX))
We call such Φ and f Φ a triangulation of the Pareto set (f) and the Pareto front fX*(f), respectively.
2.3 Bézier simplex
Let N be the set of nonnegative integers (i.e., N (0,1,2..... and
No - =
-((d.....da))) dM) ENM
Σ dm - dm= m==1
For t: = (t))...,tM) tM) A-1 and d= (d1....,dy, dM) E eN we denote by ta monomial 14 ... ...MM The Bézier simplex of degree D in RR_ with control points (pa)deNy C RL is map b: A R&, which is defined by
b(t P) = Σ
à
t Pd:
(1)
dENN
D where (2)= dy! is multinomial coefficient and P E R |o which is defined as
represents matrix of control points,
(pt)1
(P1)2
(p1)
(p:)1
(P2)2
(p2)1
P =
(2)
(P ppyB) )
(p Payl )2
(p )r
Define z(t) as vector of cocfficient of the Bézier simplex (1) with respect to control points {Pald: i.c.,
d t T
x(t) = ((0)es... tel
(p))
E R No
Then, the definition of Béxier simplex (1) is represented as = ztt.. It. is known that Béxier simplex is universal approximator of continuous functions (Kobayashi et al., 2019), and thus, the mapping argmin E(f) can be approximated by Bézier simplices arbitrary precision. From this thcoretical advantage, we construct a general framework to obtain a multi-objective oplimization method from single-objective optimization method with Bézier simplices.
PT
3 Proposed Algorithm
A number of methods have boen studied in the context of multi-objective optimization. Many of the methods are designed to apply to any multi-objective optimization problem; however, the individual methods are written in separate contexts and are not unified. Therefore, in this paper, introduce general framework to obtain a multi-objective optimization method M(A) from single-objective optimization algorithm A. Moreover, in contrast to the existing methods, which aim to find a finite set that approximates the Pareto set/front, our proposed hothod obtains a parametric hypersurface representing the Pareto eeo/front of multi-objective problems to solved.
In our proposed algorithm, we aim obtain control points of a Bézier simplex that represents the Pareto set of a multi-objective problem from single-objective optimization algorithm A. In this paper, single-objective
4

---

Published in Transactions on Machine Learning Research (04/2025)
p(tii
b(tkii Auri Arp (pii)
Aubb(ty) 1 ppk)
x
4
x2
X2
X2
X2
(a) Generate solutions the Bézier simplex.
(b) Update each solution with (0
(c) Update the Bézier simplex.
Figure 2: Conceptual diagram of M(A) the kth iteration; the red surfaces in (a) and (c) represent the Bézier simplices.
optimization algorithm A is a map from the direct product of the sample space Z and the space of loss functions C to the space of model paramcters W, i.c., A: 2 x L - W. Then, for any t E A M-1,, we denote by At an update rule in A, which defined by the loss function [x((F)
Our algorithm begins by setting the initial control points p(),, At the k-th iteration (k 2 1), we randomly draw N samples (tn J from the uniform distribution on A M-1 and obtain N data points (b(th (8) [ [
the Bézier simplex defined by the current control points p(k) Next, update each b(th.) i.e.,
(k) p(k)) by A,00,
(k)
Agg (b(t!)) p()))).
(3)
Then, we update the Bézier simplex with the obtained pairs of data {(th") whh a > Specifically, solve the following least squares problem to fit a Bézier simplex to {(th") (2) ) >
Σ 20 b(t(h) p) PER Ni XL
minimize
(4)
where P is a decision variable to be optimized, and 11.11, denotes the Euclidean norm on R& Let (+) (k)
= (o(t), 2² .......,er ( e NXL and (k)) = (x(t(()), ), (t((() ),...,, ( (T Nx|Ngl be matrix of (1) and z(thh)) respectively. Then, the problem (4) is equivalent to the following problem:
minimize Pe IL
1 N
(x) 2 200pp p
2
where I| denotes the Frobenius norm on RNLL_ Since the optimization problem (5) is an unconstrained convex quadratic optimization, and it can be shown that the symmetric matrix (k)T ((k) is regular with probability 1 (refer to Lemma B.1 and its proof), the update rule for control points described follows:
p(k+1) = (T) z(k)
Z (k)
x (4).
(6)
We repeat this procedure until k reaches the maximum number of iterations specified by the user. We summarize the pscudocode for solving multi-objective optimization problems with the above procedure in Algorithm 1 and show its conceptual diagram in Figure 2.
4 PAC Stability and Generalization Gap
Assume that there is an unknown distribution D over some space 2. We take S = (L,,..., (N) of N examples drawn iidd from D. Then, the generalization error is defined by
R
EeND (((P t)),
5

---

Published in Transactions on Machine Learning Research (04/2025)
Algorithm Multi-objectiive Optimization Method M(A) from Single Optimization Method A
Set k +-1 and the initial control point p(k)
while & S K do
3: Draw (k)
In=i for which cach (*) is drawn ii.d. from the uniform distribution on A A-11 Obtain (b(th")
Obtain ( (1) |P()))) by Equation (1).
J==1 by Equation (3).
6: Update control points by Equation (6).
7:
t k +1. 1.
end while
return p(K+1),
where E L is a given loss function, and { () denotes the loss of the model described by P with input t. Since the generalization error cannot be measured directly, we instead consider the empirical error defined by Ps = N Σ {(P) t... Then, the gencralization gap of P is defined as the difference between empirical error and generalization error, i.e.,
Rs [P] R[P].
We consider a potentially randomized algorithm A (e.g., stochastic gradient descent) and the expected value of Equation (7):
[x(Rg[A(S)] R[A(S)),
(8)
where we represent A(S) as E) for notational simplicity.
treat the approximate behavior the expected value with respect to the sample, we consider the following. First, take an event C c z that has a high probability of occurring. Then, the conditional generalization
error under the condition C' is defined by:
N
(( u.)
PPP] = B
where DN is the conditional probability distribution of C. Note that if C = NN R(P) is equal to R(P).
Next, we consider the approximate expected value of Equation (8) by
B& Rs [s())()) A(((S))
where És is the conditional expected value of C. This invariant allows discuss the expected value the generalization gap with respect to events. The following introduces the definition of probably approzimalely correct (PAC) uniform stability. This is PAC-like expansion of the uniform stability in Hardt et al. (2016).
Definition 4.1. randomized algorithm A PAC uniformly stable if for any e E (0, 1), there exists and an event De C g+1 which occurs with probability at least 1 1- such that
sup (AA E(A(S)) O) <88,
(10)
where S (li,..., IN) and S' (L,...., ..... Iv) samples differing in most example, drawn from D. satisfying (t,...., ., t, th 1.41..... IN) D.- Furthermore, PAC uniformly stable randomized algorithm A is decomposable if for any s e (0, 1), there are events B, C Z such that D. BN++
With the PAC stability, we show the following, which ensures that if an algorithm is PAC uniformly stable, the difference between its generalization and empirical error small with high probability.

---

Published in Transactions on Machine Learning Research (04/2025)
Algorithm Surface-wise Gradient Descent Method
1: Sct h <- and the initial control point P(k).
2: while & S do
3: Draw (
In=1 for which each (/ is drawn i.dd from the uniform distribution on M-1
Obtain {b(th)) (k) [P(())) by Equation (1).
5: Obtain fan In=1 by Equations (3) and (11).
Update control points by Equation (12).
k + k 1.
end while
return p(K+))
Theorem 4.2 (Proof is shown in Appendix A). Let A be a decomposable PAC uniformly stable randomized algorithem. Then, for any & E (0, 1) and 8 > 0 Theorem 4.1, there erists an event C C ZN which occurs with probability at least 1 such that,
EsEA [ALA(S]]-R[AASS <8.
where És is the conditional expected value of C. and R is the conditional generolization error under the condition Ce.
5 A Surface-wise Gradient Descent Method
Next, we discuss the case that an algorithm A is a gradient descent method. We refer the method as surfacewiise gradient descent method since the method is designed to generate sequence of hypersurfaces. Here, we employ the gradient descent-based update rule in Equation (3) as follows:
A,00 ((((()) p(k)) = ((t(()) ((l) - a (k) (b(t(h) p(k)) (())
(11)
T ((k)
b(t p(()) (*)) ((t(())
(R) p(())
where A,g) is the update rule, (*) (0,1) is step size, da is first derivative with respect to æ, f(-1t) is weighted sum of objective functions fi...., fM by t, and (y(#) is a matrix gradient fm at x defined by (x( (Vfi(z)..... ..., (fa( e MXL, Define (()) and ((k) as
B(a) (()) G( ((11)) J3 (p(P 2(14))) ..... (B J J (P(N)T((11))))
(
Note that (() and (k) are variables determined by (1) which is drawn i.idd from the uniform distribution on A M-1 in cach iteration, however, the argument abbreviated for the sake of simplicity. Then, the update rule with gradient descent method described in Equation (11) is rewritten as
kk)
(
(N) (R)
With this notation, the update rule for the control points in Equation (6) is represented
p(k+1) p(k) (k) (z Z 2(1) Z ( ((k)
(12)
We summarize the pseudocode the surface-wise gradient descent method in Algorithm 2.
6 PAC Stability of the Surface-wise Gradient Descent Method
We prove that the surfaco-wist gradient descent PAC uniformly stable. All omitted proofs are shown in the Appendix. Hereinafter, make the following assumption about the objective function.

---

Published in Transactions on Machine Learning Research (04/2025)
Assumption 6.1. All the objective functions f1...., . far are p-Lipschitz continuous and differentiable on X.
Let x": A:AA--1 + (f) be map from A M- to the Pareto set of f. For P defined in Equation (2), we define a loss function as [(P))) = [|b(t P) -)-x+"())) (t)l Since (f) is unknown, we cannot take sample directly from (f). Instead, we take sample [t- drawn i.i.d. from the uniform distribution on A
To prove that Algorithm 2 is PAC uniformly stable, we show two propositions in advance. Note that the following two propositions can respectively be regarded as an extension of boundedness and expansiveness introduced in (Hardt et al., 2016) to analyze the stability of an optimization algorithm.
Lemma 6.2. Let U > O be constant satisfying maxtes"-11 ||((()]]]]] SU Let Pr be the update rule with parameter T = [tn]] j in Equation (12), i.e., p(++1) = (r (P())) Then there exists some n > O, and we have the following inequality with probability at least 1 -: E:
1 (P) - Pllp S NNUU..
Lemma 6.3. Let > O and U > 0 be constants as in Lemma 6.2. For T = (t and T {t Ji such that the difference between T and T lies only in one example, there crists some Ç > 0, and we have the Sollowing with probabilily al least 1 -
(pr(())rr(PPll) S pU (N).
Let (T: and {T be parameters whose difference lies only in the k-th element, and p(K+1) and p((K+l) be respectively the output of Algorithm 2 with {Ti) and (T/) From Lemmas 6.2 and 6.3, we show that |p(k+1 - p* is bounded above with arbitrary probability.
Lemma 6.4. Let n and <> be constants as in Lemma 6.3. Suppose that run Algorithm 2 for K iterations with paramelers {T: and (T whose difference lies only in the kth element. Then, we have the following with probability at least 5:
2qnn (+ [ Ck-k- k+ ))x))
- - picktt)|] -piKktll S
Now, we are ready show that Algorithm 2 is PAC uniformly stable. Theorem 6.5. Assume that (k) E (0,1) 1] for all [K]. Then, Algorithm 2 is PAC uniformly slable.
From Theorems 4.2 and 6.5, if Algorithm 2 is decomposable, we can obtain an upper bound of its generalization gap.
Remark 6.6. While we assume the decomposability of Algorithm 2 in Theorem 6.5, this assumption is not restrictive in practice for the following reason. Suppose that an Algorithm 2 is not decomposable, meaning that there exists ह E (0, 1) such that there are no event BN satisfying D. BB Even in this case, we can often take a smaller z' < so that there exists an event Be satisfying B B C Dz. Thus, by letting an event as De + BN++1 we can bound the difference between the generalization and empirical errors with an event De which occurs with a certain probability.
7 Numerical Experiments
To verify that the Pareto optimal set can be accurately approximated by a Bézier simplex obtained by the proposed method (Algorithm 2), we applied Algorithm 2 to three multi-objective problems: scaled-MED, skew-M MEED, and MMDW employed in (Hamada et al., 2011), which are known to be simplicial. Skew-MMMD includes some important real-world problems such as group LASSO in sparse modeling (Yuan & Lin, 2006), and MMED includes generalized location problem (Kuhn, 1967). For real-world instances, we applied Algorithm 2 to a group LASSO problem with the Birthwt dataset from the MASS package" in language. In addition, to confirm the applicability of the proposed method to non-simplicial
'https://cran.rr w-wonnertrroottcoo.e. progjjtt.ogg/ee///cccnggss//eeee
html
8

---

Published in Transactions on Machine Learning Research (04/2025)
Table 1: MSE (avg.+s.d. over 100 trials) for scaled-MED.
N
Proposed (Algorithm 2)
NSGAIT + all-at-once MOEA/D /D + all-at-once
30
6.07e-05 + 3.18c-06
9.31c-02 + 5.40c-04
3.32c-01 + 4.28c-03
50
4.63e-05 + 2.41e-06
1.37e-01 + 8.18c-04
1.26c-01 + 6.80e-04
100
3.83e-05 + 1.59e-06
8.86e-02 + 7.16e-04
1.14e-01 + 6.62e-04
instances, we applied Algorithm to DTLZZ instances (Deb et al., 2002); which are widely used as benchmarks for multi-objective optimization algorithms. The definition of each problem instance is shown in Appendix E. In Algorithm 2, we set the degree of Bézier simplex as D 3, the initial control points as p()) = 0, which is the zero matrix of appropriate size. In addition, set the maximum number of iterations K 100 and the step size as (*) = 1 for k = .... K. iteration was tested for N E {30,50, 50, 1003.
The number of points to be sampled from a simplex in each
As a baseline, we used NSGA--I (Deb et al., 2000) and MOEA/D (Zhang & Li, 2007) with the Bézier simplex fitting (Kobayashi et al., 2019). Specifically, we obtained approximated Parcto solution samples by NSGA-II and MOEA/D implemented in jMetal (Benitez-Hidalgo et al., 2019), with default parameters except for population size. Then, we fitted the approximated Pareto solution samples with Bézier simplex of degree D = 3 by the all-at-once method proposed in (Kobayashi et al., 2019). We set the population size N E {30,50) 1003. We implemented these algorithms in Python 3.12.3, and the experiments were performed on a Windows 10 PC with an Intel(R) Xeon(R) W-1270 CPU 3.40 GHz and 64 GB RAM. The source code is accessible at https: //github. www.c/oooooooi//eeeeeeeeeeeaaaee hiw.inaa/bbziie--fiioo
7.1 MSEs comparison
First, we picked up simplicial problem instance whose map xt: A-1 - X*f) is analytically obtained and evaluated how accurately an obtained Bézier simplex approximates the optimal Pareto set. In this experiment, we used scaled-MED, which is a three-objective problem with three variables and is known to be simplicial. The problem definition is shown in Appendix E. To evaluate the approximation accuracy of the estimated Bézier simplex, we used the mean squared error (MSE) defined by MSE = &/ El=,llb(En P) - x² (t..) l| where x" is a map from a weight t E A to the minimizer (t) of the corresponding scalarizing function. The map x² for scalod-MED is shown in Appondix F. To calculate MSE, we randomly sample {f, ) from the uniform distribution on 2 We repeated the experiments 100 times with different parameters and computed the average and the standard deviations of MSFs.
10000 i.i.d.
Table shows the average and the standard deviation of the MSES with N (30,50, 50, 100) for Algorithm 2 and population size N E {30, 50, 100} for NSGA-II and MOEA/D. In Table 1, we highlighted the best score of MSF out of the proposed and baseline methods where the difference is significant with the significance level 0.001 by the Wilcoxon rank-sum test. Table 1 implies that the Bézier simplex obtained by our proposed method can represent the Pareto set well. Moreover, the MSEs of our mothod decrease with larger N, which supports the PAC uniform stability of Algorithm 2.
Figure 3 shows the Bézier simplex obtained by our proposed method, and Figure shows the Bézier simplex obtained by the all-at-once with the approximated Pareto solutions of NSGA-II. The true Pareto set of scaled-MED is known to be a curved triangle that can be triangulated into three vertices. Recall that the analytical solution of scaled-MED is shown in Appendix F. In Figure 3, the Bézier simplex obtained by the proposed method approximates the Parelo set well even when N = 30, while the Bézier simplex obtained by the all-at-once method with NSGA-II does not approximate the Pareto set even when N = 100. As a supplementary empirical ovaluation, we conducted experiments to investigate the algorithmic stability of the proposed method on the scaled-MEN The details and results of this experiment are presented in Appendix H.
9

---

Published in Transactions on Machine Learning Research (04/2025)
Bézier simplex
Optimal Pareto set
Approximated Pareto set
10
1.0
1.0
0.5
0.5
0.5
0.0
0.5
005
-0.5
-100
-1.0
-1.0
0.90.20.4 A..60.8
600.0.60.40.0
0.9..7 0.6
600.00.60.0022
0.00. 0.00.70.4
1.0.00.604002
20.00
0.60.81
22
21
11
23
3
5
-
Figure 3: Results for Algorithm with the sample size N E {30, 50, 100}.
Bézier simplex
Optimal Pareto set
Approximated Pareto set
1.0
1.0
1.0
0.5
0.5
0.5
0.0
0.0
0.0
0.5
-05
0.5
-1.0
0.2
0.6 0.4 0.2000
0.0
0.2
0.2
0.6 0.4 0.2
42
0.4
0.6
0.8
1.0 1.0
42
0.4
0.6
1.0
1.0
0.8
0.6
0.4
21
32
0.4
0.6
0.8
1.0
1.0
Tx
N
30
N 50
N 100
Figure& 4: Results for NSGA-II and all-at-once with the population size N {30, 50, 100}.
7.2 GDs and IGDs comparison
Next, we validate the practicality of the proposed mothod in more practical scttings. In this experiment, we used two synthetic simplicial problem instances (skew-3MED and skew-3MMD) and one real-world problem instance (group LASSO), whose map x* cannot be represented in closed-form. Additionally, to confirm the applicability of the proposed method to non-simplicial instances, we applied Algorithm 2 to DTLZ instances (Deb et al., 2002) as non-simplicial instances. We show the definition of each problem instance in Appendix E. We used the generational distance (GD) (Veldhuizen, 1999) and inverted generational distance (IGD) (Zitzler et al., 2003) to evaluate how accurately the estimated Bézier simplex approximates the Pareto optimal set, which is defined as follows:
GD(X,Y) =
min | vlla. IGD(X,Y) Σ mEX
Σ BEY
min - y12.
X
yeY
Y
where X is finite set whose elements are sampled from an estimated hypersurface and Y is a validation set. We can say that the obtained Bézier simplex is close to the Pareto set if and only if both GD and IGD are small. As validation set Y, generated approximate Pareto solutions by NSGA-III (Deb & Jain, 2013; Jain & Deb, 2013) with the population size of 1000. To construct X, we randomly sampled 11000 i.i.d. from the uniform distribution on A and obtained sample points on the estimated Bézier simplex. We repeated the experiments 100 times with different, parameters and computed the average and the standard deviations their GDs and IGDs.
Tables 2 show the average and the standard deviation of the GDs and IGDs for 3kew-3MED, Skew-3MMD.. and group LASSO with the number of samples N (30, 50, 100) and the population sizes N (30, 50, 100).
10

---

Published in Transactions on Machine Learning Research (04/2025)
The results for non-simplicial instances are shown in Appendix G. In Tables 2 to 4, we highlighted the best score of GD and TGD where the difference is at significant with significance level p = 0.001 by the Wilcoxon rank-sum test. Tables 2 to 4 show that the proposed method achieved better GD and IGD for all cases. The differences are pronounced in the results of small sample/ sample/population population size, which implics our method obtains a Bézier simplex approximating Parelo set well.
7.3 Run time comparison
To evaluate the computational elficiency, we measured the computation time of the proposed method (Algorithm 2) and two baseline methods (NSGATT and MOFA/D. with all-at-once) for skew-3MED and skew-3MMD. Table 5, we show the average and the standard deviation of the computation time in seconds of each method over ten trials. For the proposed method, we report the computation time with the maximum number of iterations K = 100. From Table 5, we see that the proposed method took about ten seconds to output the resulting Bezier simplex, whereas the post-processing method (NSGA-1L and MOEA/D with Bezier simplex fitting) needed only about one to two seconds for each N. This is mainly because our method performs Bezier simplex fitting at every iteration, whereas existing methods do this only once. On the other hand, the computation time for the proposed methods remained almost the same even when N increases by 100, while it grows for the existing methods. Although our method was slower in our experiments, its computing time did not increase much with larger N, which can be advantageous in certain situations.
8 Conclusion
This paper proposes a general framework to construct a multi-objective optimization algorithm from single objective method with the Bézier simplex. We have also defined the PAC stability of optimization algorithms and proved that this stability gives us an upper bound on the gencralization gap in the sense of PAC learning. The theoretical analysis showed that if we construct multi-objective optimization algorithm from a gradient descent-based single-objective optimization algorithm, the resultant algorithm is PAC stable. Tn numerical experiments, applied the multi-objective optimization algorithm with our framework to three synthetic and one real-world instances and demonstrated that our algorithm attained better gencralization gaps and approximation ion accuracies of the Pareto optimal set than the existing algorithm. As a concluding remark, we have to note that this study is limited to treat simplicial problems. It would be interesting for future studies to extend this study to non-simplicial cases.
References
Antonio Benítoz-Hidalgo, Antonio J Nebro, José García-Nicto, Izaskun Oregi, and Javier Del Scr. jmetalpy: A python framework for multi-objective optimization with metaheuristics. Swarm and Evolutionary Computation, 51:100598, 2019.
Sergio F. Contreras, Camilo A. Cortés, and María A. Guzmán. Bio-inspired multi-objective optimization design of highly efficient squirrel cage induction motor. In Proceedings of the Genetic and Evolutionary Computation Conference 2016, pp. 549 556, 2016. doi: (0.1145/2908812.290
Kalyanmoy Deb. Multt-objective Optimization Using Evolutionary Algorithms. John Wiley & Sons, Inc., New York, NY, USA, 2001.
Kalyanmoy Deb and Himanshu Jain. An evolutionary many-objective optimization algorithm using reference- point-based nondominated sorting approach, Part Solving problems with box constraints. IEEE Transactions on Evolutionary Computation, 184))5577-001 2013.
Kalyanmoy Deb, Samir Agrawal, Amrit Pratap, and Meyarivan. A fast clitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II In International Conference on Parallel Problem Solving from Nature, pp. 849 858-858.866 Springer, 2000.
11

---

Published in Transactions on Machine Learning Research (04/2025)
Table 2: GD and IGD ts.g. over 100 trials) for skew-3MED.
N
Proposed (Algorithm 2)
NSGAIT + all-at-once
MOFA/D + all-at-once
GD
30
7.07c-02 + 2.05c-03
2.59c-01 + 7.41c-03
1.57c-01 + 5.42c-03
50
6.53e-02 + 1.18c-03
1.26c-01 + 4.72c-03
9.36c-02 + 1.46e-03
100
6.59e-02 + 1.23e-04
7.57e-02 + 1.11e-03
6.63e-02 + 9.19e-04
IGD
30
9.91c-02 + 8.43c-04
1.69c-01 + 5.39c-03
1.28c-01 + 3.35c-03
50
9.75e-02 + 4.84e-04
1.15e-01 + 1.89e-02
1.09e-01 + 3.47e-03
100
8.98e-02 + 5.22e-04
9.19e-02 + 3.44e-03
9.77e-02 + 2.78e-03
Table 3: GD and IGD (avg. ts.d. over 100 trials) for skew-3MMD.
Proposed (Algorithm 2)
NSGAII + all-at-once
MOEA/D + all-at-once
GD
30
5.57e-02 + 1.69e-03
1.61e-01 + 3.19e-03
5.94e-02 + 9.98e-04
50
5.33e-02 + 1.25c-03
9.98e-02 + 1.98c-03
5.94c-02 + 9.86e-04
100
5.14e-02 + 1.02e-03
8.57e-02 + 2.50e-02
5.94e-02 + 2.54e-03
IGD
6.69e-02 6.73e-04
8.57e-02 +: 2.78e-03
8.58e-02 +: 3. 3.16e-03
50
6.34e-02 5.57e-04
7.17e-02 + 2.48e-03
6.45e-02 = 3.06e-03
100
6.28c-02 + 4.62o-04
1.01c-01 + 6.31c-02
7.31o-02 + 2.95c-03
Table 4: GD and IGD is.d. over 100 trials) for group LASSO.
N
Proposed (Algorithm 2)
NSGAII + all-at-once
MOEA/D + all-at-once
GD
30
2.44e-02 + 5.25e-04
6.47e-02 + 1.75e-04
6.78e-02 + 2.78e-03
2.37e-02 4.72e-04
4.29e-02 +: 3.73e-04
4.13e-02 +: 4.51e-04
100
2.37e-02 t 4.30e-04
4.56e-02 + 5.06e-04
4.12e-02 + 3.10e-04
IGD
30
5.94e-02 + 3.98e-04
2.05e-01 + 4.12e-04
1.35e-01 + 2.76e-04
50
5.95c-02 + 3.79c-04
1.18c-01 + 3.99e-04
9.96c-02 = 3.97e-04
100
5.92e-02 + 3.82e-04
1.24e-01 + 5.13e-04
6.44e-02 + 3.12e-04
Table 5: CPU time in seconds for skew-3MED and skew-3MMD.
N
Proposed (Algorithm 2)
NSGAII + all-at-once
MOEA/D + all-at-once
skew-3MED
30
1.09e+01 5.79e-01
1.41e-00 + 1.09e-01
1.92e+000 4.92e-02
1.07e+01 + 7.25e-01
1.67e-00 4.0Ge-02
2.04e+000
5.38e-02
100
1.08e+01 = 6.72e-01
2.57e-00
8.76e-02
1.85e+00
3.74e-02
skew-3MMD
30
1.32e+01 + 7.26e-01
1.39e-00 + 8.89e-03
2.04e+000 +:2 5.38e-02
50
1.29e+01 = 5.84e-01
1.83e-00 = 9.01e-02
2.12e+000
2.34e-02
100
1.30c+01 6.95c-01
2.69c-00
7.59c-02
2.43e+000 2.89c-02
Kalyanmoy Deb, Samir Agrawal, Amrit Pratap, and T. Meyarivan. A fast and elitist multiobjective genetic algorithm: NSGA-IL. IEEE Transactions Evolutionary Computation, 6(2):182 197, April 2002. doi: 10.1109/4235.9900
Gabriele Eichfelder. Adaptive Scalarization Methods in Multiobjeclivee Oplimization. Springer-Verlag. Berlin, Heidelberg, 2008.
Naoki Hamada and Shunsuke Ichiki. Simpliciality of strongly convex problems. Journal of the Mathemationl Society of Japan, 73(3), 2020. doi: 10.2969//mmsj/8331888
Naoki Hamada, Yuichi Nagata, Shigenobu Kobayashi, and Isao Ono. On scalability of adaptive weighted aggregation for multiobjective function optimization. In 2011 IEEE Congress of Evolutionary Computation (CEC), pp. 669 678, 2011. doi: 10.1109/CEC.201 CFC.2011.599955
12

---

Published in Transactions on Machine Learning Research (04/2025)
Naoki Hamada, Kenta Hayano, Shunsuke Ichiki, Yutaro Kabata, and Hiroshi Teramoto. Topology of Pareto sets of strongly convex problems. STAM Journal on Oplimizalion, 30(3):2659-2688 2020. doi: 10.1137/19M1271499.
Ken Harada, Jun Sakuma, Shigenobu Kobayashi, and Isao Ono. Uniform sampling of local Pareto-optimal solution curves by Parcto path following and its applications in multi-objective GA. In Proceedings of the Genetic and Evolutionary Compulation Conference 2007, pp. 813 820. ACM, 2007. doi: 10.1145/1276958. 1277120.
Moritz Hardt, Ben Recht, and Yoram Singer. Train laster, generalize better: Stability of stochastic gradient descent. In Proceedings of the 33rd International Conference on Mochine Learning, volume 48, pp. PMLR, 2016.
Daniel Hernandez-LLoatto Jose Hernandez-Lobato, Amar Shah, and Ryan Adams. Predictive entropy search for multi-objective bayesian optimization. In Proceedings of the 33rd International Conference on Machine Learning, volume 48, pp. 1492-1501. PMLR, 2016.
Claus Hillermeier. Nonlinear Multiobjective Optimization: A Generalized Homotopy Approach, volume 25 of International Series of Numerical Mathematics. Birkhäuser Verlag, Basel, Boston, Berlin, 2001. doi: 10.1007 97000099998---34-828---2
Himanshu Jain and Kalyanmoy Deb. An evoautionary many-objective optimization algorithm using reference- point based nondominated sorting approach, Part II: Handling constraints and extending to an adaptive approach. IEEE Transactions on Evolutionary Computation, 18(4):602 622,,26222222222 2013.
Ken Kobayashi, Naoki Hamada, Akiyoshi Sannai, Akinori Tanaka, Kenichi Bannai, and Masashi Sugiyama. Bézier simplex fitting: Describing Pareto fronts of simplicial problems with small samples in multi-objective optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pp. 2304-23112 2019. doi: 10. 1609/ (6609/aaaivva. a3ai.v33:01.3301223
Harold W. Kuhn, On a pair of dual nonlinear programs, Nonlinear Programming,
1967.
Stefanus C. Maree, Tanja Alderliesten, and Peter A. N. Bosman. Ensuring smoothly navigable approximation sets by bézier curve parameterizations in evolutionary bi-objective optimization. In Parallel Problem Solving from Nature PPSN XVI, pp. 215 228. 2020. URL
Franco Mastroddi and Stefania Gemma. Analysis of Pareto frontiers for multidisciplinary design optimization of aircraft. Aerospace Science and Technology, 28(1):40-55, 2013.
doi: 10. 10.1016/jaas.20122.0000
Kaisa Miettinen. Nonlinear Multiobjective Optimization. Springer US, 1999. doi:
Yusuke Mizola, Naoki Hamada, and Shunsuke Ichiki. All unconstrained strongly convex problems are weakly simplicial. arXiv, arXiv:2106.12704, 2021.
Gilberto Reynoso-Meza, Leandro dos Santos Coelho, and Roberto Z. Freite. Efficient sampling of pi controllers evolutionary multiobjective optimization. In Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation, pp. 1263-1270, New York, NY, USA, 2015. Association for Computing Machinery. ISBN 9781450334723.. doi:
Siarnak Safarzadegan Gilan, Naman Goyal, and Bistra Dilkina. Active learning in multi-objective evolutionary algorithms for sustainable building design. In Proceedings of the Genetic and Evolutionary Computation Conference 2016, pp. 589-596, 2016. doi: 10.1145/2908812
Oren Shoval, Hila. Sheftel, Guy Shinar, Yuval Hart, Omer Ramote, Avi Mayo, Erez Dekel, Kavanagh Kavanagh, and Uri Alon. Evolutionary trade-offs, Pareto optimality, and the geometry of phenotype space. Science, 336(6055):1157 1160, 2012. doi: 10. 10.1126/science. 1217400..
Akinori Tanaka, Akiyoshi Sannai, Ken Kobayashi, and Naoki Hamada. Asymptotic risk of bézier simplex fitting. In Proceedings of the AAAI Conference on Artificial Intelligence, pp. 2416-2424 2020. doi: 10. N0.1609/aaii.v34103..66
13

---

Published in Transactions on Machine Learning Research (04/2025)
David Allen Van Veldhuizen. Multiobjective Evolutionary Algorithms: Classifications, Analyses, and New Innovalions. PhD thesis, Department of Electrical and Computer Engineering, Graduate School of Engineering, Air Force Institute of Technology, 1999.
Jasper A. Vrugt, Hoshin V. Gupta, Luis A. Bastidas, Willem Bouten, and Soroosh Sorooshian. Effective and efficient algorithm for multiobjective optimization of hydrologic models. Water Resources Research, 39(8): 1214-1232, 2003. doi: 10. 20.1029/200WRR001722
Yieh-Hei Wan. On the algebraic criteria for local Pareto optima- L. Topology, 16:113 117, 1977. doi:
Yich-Hei Wan. On the algebraic criteria for local Parcto optima. II. Transactions of the American Mathematical Sociely, 245:385 397, 1978.
Kaifeng Yang, Michael Emmerich, André Deutz, and Thomas Bäck. Multi-objective Bayesian global optimization using expected hypervolume improvement gradient. Swarm and Feolutionary Compulation, 44:945 956, 2019. doi: 10.1016/j.swevo.2011 20w.::............0...0
Ming Yuan and Yi Lin. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(1):49-67, 67, 2006.
Qingfu Zhang and Hui Li. MOEA/D: multiobjectivee evolutionary algorithm based on decomposition. IEEE Transactions Fvolutionary Compulation, 11(6):712-783 731, 2007. doi: (0.109/TEVV.2200..52777
Eckart Zitzler, Lothar Thiele, Marco Laumanns, Carlos M. Fonseca, and Vviane Grunert Da Fonseca.
Performance asscssment of multiobjective optimizers: An analysis and review. IEEE Transactions on Evolulionary Compulation, 7(2):117-133) 2003. doi: 10.1099/TEVC.2003..111
A Proof of Theorem 4.2
Proof- We denote two independent random samples by S (t1 (f1,---, tN), S' (1/,...,, (ty) Let 5( (t1...., .....,t-1,66 E2 fi+1...... tN) be the sample that is same as S except in the ith example where replace t; with t.. Since D. BN+1, (t)...., t, t. t...t......... (x) E De for any i if and only if S, S' E C. B, In
this case, we have
E ((A(S) -((()))))) <
Then, adding the incqualities for and applying the triangle incquality, we obtain
N Σ ((A(S() (.))
BA
((A(S)
<8.
N
Let us denote the conditional probability distribution of D, DN under the condition B2, C. by B Ce respec- tively. Then we have
N ((A((())))))
((A(S)
N Σ
<6.
N
1=1
Here, we have
E E( EA
((A(S) (c:)
Es-cEE.
Σ ((A(S)
N
N
1=1
14

---

Published in Transactions on Machine Learning Research (04/2025)
and
N
NV Σ
C2
Σ ((A(S()) [t,) 2m1
((8.8))~c3BA CA
((A(S) t
N
.V
N
Es BS
Σ 1=1
((A(S))E)
N
EsEA R[A(S)]]
where Es is the conditional expected value of C. and R[A(S) is the conditional generalization error of C..
Thus, we obtain the inequality in the theorem.
Finally, we have
P(C,) = P(B)) P(D:) Ai P(Ds)
which completes the proof.
B Proofs of Lemmas 6.2 and 6.3
We first show the three lemmas in advance. Lemma B.l. For all (), there crists v) > O satisfying
P Pmij Amiin (TT)) Z(T.)) 21
where Amin denotes the minimal eigenvalue of symmetric matria A, and T; = (thh)) LN from the uniform distribution on A M-1 for [K + 1].
Jn=1 is drown i.i.d.
Proof. Consider the set
Fo = [ S (AM-11)
K+1
min Amin (Z(T)) =0)
Since Z(T.) (TT)) is a symmetric matrix, has non-negative eigenvalues and Z(T)) Z(T)) has zero cigenvalue if and only if the determinant of Z(T:) Z(T:) is zero. Hence, Fo is subsct of
(tn) C-M-1 M-1 Il del Z(T:) Z(T)) 0
Therefore, Fo is equal to the zero set of polynomial. This implics P(Fo) 0. Considering the set K+1 min Amiu (T)))) T z(T.)) S $n}
Fy = (T. C (AM-1)K*
then we have
P(F,) - P(Fo) = 0 (n + 0).
This implics that for all E (0, there exists some 1 > 0 such that
P min Amin (TT)) (TT)) sn) E.
(13)
The proof is completed by taking the complementary event.
Lemma B.2. For all e E (0, 1), there erists Ç > 0 satisfying
IP
max
( Z(T))
<c)
(14)
where T. ( Jn=1 is drawn i. from the uniform distribution on A-1 for E [K 1).
15

---

Published in Transactions on Machine Learning Research (04/2025)
Proof. Since (T;) ((T)) is symmetric matrix, there exists orthogonal matrix Q: such that (T:) ( = Q where A, is a diagonal matrix whose diagonal entry is the eigenvalue of ( (T)). According to Lemmn B.1, (TT)) Z(T)) is a regular matrix with probability at least 1 - e for all i e [K + 1]. Hence, we have the following with probability at least -:
max (Z(T) ((T))
max Q A Qullr
max
AA7 I|
max
X Z(T,) Z(Ti)
ND
min;
X2 min
( T (
NW
where the second equality follows from the fact that the Frobenius norm is unitarily invariant, and the second incquality follows from Lemma B.1. The proof is completed by letting be some real number greater than or equal to V
1
Lemma B.3. Let > 0 be a constant satisfying =axxesM==11=(( U. Then, for any S M-1 have
we
| T G
S NUn.
(15)
Proof. First, we show that ||=())]2 is bounded above for any t AM Since is continuous function
over t whose domain A M-1 is compact, there exists upper bound U > O for any t E A A-1, Next, we show that F
is bounded above. For any T = (t C A M M-1, we have
JP 21)
(TT) T G(T) = (211..., .N)
zn t (PT
2x)
n=1
Therefore,
N
Z(T) T G(T)
Σ
nn 1 J (P (PTzz)
=
n=1
P
2nt, Jj (PT 2n)]l
< >
|| 2 1
|| Jy 2,1
N
M
Σ
|| ||
Σ m=l
Imm V (PT
n=1
<
Up NUp.
n=1
16

---

Published in Transactions on Machine Learning Research (04/2025)
The last inequality holds by the fact that the term Em= Fm Fm (P z)...... ...., Vf P z) and the assumption that every function fi,..., s.....faa fM u-Lipschitz continuous.
(P 2) is convex combination of
Finally, we show Lemmas 6.2 and 6.3.
Proof of Lemma 6.2. By Equation (12), we have
11pr(P) - pllp = (bk)
(Z (z 2) -1 (3 ZTG
F
S (k)
(z z) 1|
z G|ll l
Let > 0 be a constant as in Lemma B.2. From Lemmas B.2 and B.3, we have the following with probability
at least 1 - s:
[ppe(P) -Plpp S NUU,.
Proof of Lemma 6.3. Let T, T" and T be
T {11 ....... 1
T 141 1N L' J
T (t.....,ty.. IN-I, t'y
Let Z be matrix constructed by T. By Sherman-Morrison formulaa, we have
( 2 (z Z #NN-1
1 Z) ZN+12 N+1 (Z (ZTZ)" 2)
=z=
1 N 2) 2N+1
Let up (P) be the control points obtained by Algorithm with T. Then, we have
(P((P) -p)() (
Z - (z ZG
22 1 - - zT 2z
2N 2 N+1 (ZT (ZTZ)= - ZTG
1 1+2 [2 (ZTZ) 2N+1
(z Z) 2N+1*N+11 (Z 2) ZT
(Z -Z ( 2N41))
1 x x (Z) 2) 2N+1
Considering the norm on both sides, we have
(P) (r(()) s <|(
|| || (PT z ll2
(2TZ) T 2 (2
+ Z G
F
1+ xN41 (Z TZ)) 2N+1
In the following, for the sake of simplicity, let A = T Z, = zN+1 and = A-1b. 16. Then, we have the following incquality with probability at least - 5:
2 2N+12 1 2)
A A- A
+ TT [ 2TZ)" 2N+1
FF
1+ bT A b
F
17

---

Published in Transactions on Machine Learning Research (04/2025)
( A (b A A 1+b A-1 4AAAAA b
yy
1 + B Ay |F
1 | 1
Ayy
1 12
y Ay
Amin (A) < 5:
where Ç is constant in Lemma B.1. The first inequality holds since := ZT Z is a positive semidefnite matrix with probability 1 by Lemma B.1, and the second inequality follows from the property of Rayleigh quotient. The last incquality directly follows from Lemma B.1. Hence, we have the following inequalities with probability at least 1 8:
p7((P) pT(P) S <1 (N),
and
(pr(P)))r(PP)) (r)l S (n
Therefore, we have the following with probability at least -
[1oT(P) pr P)]l= = | PP) 1p7 + (P((P 1o (P)
$ (p7(P) -pr(P)) + +1|07(P) -p (pr (P)
< 21) + CN).
C Proof of Lemma 6.4
Proof. Let 80) p(i)_ p(()) Ipp. We have 800) for i I...., k. From Lemma 6.2, have the following with probability at least 1 -
(++1)
p((+1) p((+1)
p((+)) p(() p(t) p((t) p(()) p((L-1)
p(+1)) p(il F + pr(it1) -1041)_prhh - + +||p(0
pro
- Nn 8(1)
for each i = K...,,. Therefore, by using the above relation repeatedly and from Lemma 6.3, we have the following with probability at least 1 E:
g(K+1) 2 k)NNUp pu +NN
2 (1+ K 1K+ (k)
which completes the proof.
18

---

Published in Transactions on Machine Learning Research (04/2025)
D Proof of Theorem 6.5
Proof. For any t E A M-1 and for any (Ti) (T C (AM-1) such that (Ti) and (T differs only one example, we have
(((TT)) ((A(T)) -(A(()))))((((T)
b(t p(K+1)) x² (t)
/((t) p((k+t)) x x=((|1
bit p(k+1)) b(t) p(( 1
(16)
p(k+1) p)(ktt)) T x(t)
- [z((1
pik+n -
where the first inequality follows from the reverse triangle inequality. We can bound the right-hand hand side of Equation (16) with probability at least 1 - e by Lemma 6.4. Since the left-hand side of Equation (16) is bounded for all t E A-1, we see that Algorithm 2 satisfies PAC uniform stability.
E
Problem Definition
Scaled-MED is a vhreelvariable three-objective problem defined by:
minimize
f(x) = ((()
subject to
x E R
where
fi(x)) z1 + 3 - + 1)3,
f2(2) 2(x1 - 11)² +3
1)
fs(x) 3(x1 - 1) + 2(02 1)2 (wa 1)3.
Skew-M/ED is an M-variable M-objective problem defined by:
minimize f(x) (fi(x),....... ...., fM (m))
subject to x E RM
1
Pm
where fm(())
V2 ||00 emll
22mm - 1)
Pm exp
M
em (0,... 0, 0, .,..,,0)
m-th
for m 1,...,M. M.
Skww-MMMD an M-variable M-objoctive problem defined by:
minimize f(x) = (fi(z), far(x))) subject to where
X C R
fm (2) A (x cm
> 0,
for m = M.....MMM
In the experiments Section 7.2, we set M 3, X TR-3
3 4 4
(4 3
4
A == diag
5 Ag = diag
Ag == diag
5 5'5' 5
19

---

Published in Transactions on Machine Learning Research (04/2025)
Cm = em and Pm = exp M-1 - 1).
2(m-1) Note that diag (+) denotes the diagonal matrix.
Regarding the group Lasso instance, we followed the problem definition in (Tanaka et al., 2020, Ap- pendix E.2.2).
F Analytical Solution of Scaled-MED
We derive map x²: : A (f) for scaled-MED. For any t (2)) A?, the scalarizing function weighted by & is defined by
f(x == Σ
Im fmfm(2) m==
tΙ + 212(21 - 1)² + 1)2
+ 3/ 1)² +2 2 1)²
+ 21 (213 1)2 3/2(13 1)2 + la + 1)
Since f(=lt) is a convex quadratic function with respect to each 21, 22 and æ3, its optimal solution 2z))) satisfies the following conditions:
(f [t)
2tzz1 + 4f2 - 1) + 61z(11 1) 0,
0z
(t)
(f [t)
6t 1) + 212002 + 41/(12 - 1) = O,
Oz
(t)=((t)
05 [t)
4th (I3 1) 6tz(23 - 1) + 26z(13 1) 0.
Dzs
(t)
By solving the above equation, the map xx (t) is given by
21g + 343
3/1 + 213
21 + 312 1/
x = (z))) f1 + 2tz + 3t3" 31, + t2
2t3"
2t, + 3f2 + tz
G Numerical Experiments on Non-simplicial Instances
We present the results of Algorithm 2 for DTLZ instances in Tables 6 and 7. As in Soction 7.2, we show the average and the standard deviation of the GDs and IGDs. In this experiment, we picked DTLZ2 and DTLZ4 as non-simplicial instances. The definition of each instance can found in (Deb et al., 2002, Chapter 6.7). As can be seen from Tables 6 and 7, the proposed method does not work effectively for non-simplicial problems. Thus, the development of optimization methods for non-simplicial problems remains an important direction for future research.
H Experiments Stability
We conducted additional numerical experiments assess the algorithmic stability. Since directly evaluating
stability property is difficult, we examined the empirical difference p(K+1) - pi(K+1)|| [| in Lemma 6.4. To this aim, we repeated the following procedure ten times on the scaled-MED instance:
Generate two sets of parameters Tk = { I (1 IND...... (1 (KK)) ]i=1] and Tk =
,...,, I (E(
k = 1,2 K) is drawn 1 and (t. (k) i= are different only element. In other words, there exist [N] such that ((0) # #t,h t1000 and (()) + (() for any i j.
(K) p where each (1) (i =
|=
uniformly random from the standard simplex A M-1 Here, (t. (N)
20

---

Published in Transactions on Machine Learning Research (04/2025)
Table 6: GD and TGD (avg.ts.d. over 100 trials) for DTL.22.
N
Proposed (Algorithm 2)
NSGAIT + all-at-once
MOEA/D + all-at-once
GD
30
2.42c+01 + 1.05c-000
5.89c-02 + 1.93c-03
3.66c-02 + 4.96c-04
50
1.82c+01 = 3.69c-01
2.86c-01 + 5.28c-03
3.44c-02 = 4.35e-04
100
1.47e+01 + 2.08e-01
4.56e-02 + 3.26e-04
2.98e-02 + 3.56e-04
IGD
30
1.67c+00 = 6.27c-02
2.51c-01 + 5.59c-04
3.83e-01 + 1.83c-04
50
3.20e+000 = 2.28e-01
4.93e-01 + 6.77e-05
4.42e-01 + 1.12e-03
100
1.54e+00 = 2.64e-01
5.89e-01 + 1.40e-03
2.85e-01 = 5.71e-04
Table 7: GD and IGD (avg. over 100 trials) for DTLZ4.
Proposed (Algorithm 2)
NSGAIl + all-at-once
MOEA/D + all-at-once
GD
30
4.34e-02 + 1.56e-03
1.93e+00 + 4.39e-01
6.83e-02 + 1.75e-03
50
4.67e-02 + 1.70c-03
3.54c++00 + 9.42o-02
3.68e-01 + 5.97e-03
100
5.05e-02 +: 1.86e-03
4.55e+02 + 1.88e+01
4.07e-01 + 5.90e-03
IGD
1.97e+-00 + 2.13e-04
8.52e-02 +: 2.80e-02
4.07e-02 2.69e-04
50
1.97e+000 = 5.99e-04
6.14e-02 + 9.07e-03
4.76e-02 = 5.80e-03
100
1.97c+00 = 5.72o-04
6.51c-02 + 2.75c-02
5.41c-02 = 8.32c-03
Table 8: Empirical difference
p(K+)) p((K+])||
1p for & {20, 40, 60, 80).
== 20
=
60
= 80
p((+1) p((K+)|lpp 1.21e-03 + 5.13e-04
1.78e-03
1.10c-03
1.09c-03
+ 1.50e-03
+ 4.89e-04
+ 3.16e-04
Run the proposed method (Algorithm 2) with a maximum number of iterations K 100 for Tx and Tk and obtain the resulting control points matrix p(K+1) and p((K+)), respectively.
Calculate the difference p(K+1) p((K+))|lF.
We show the results of the average and standard deviation of | p'(K+)) for i {20, 40, 60. 80} below. From this table, we that the average of the empirical difference did not necessarily decrease for sufficiently small
large k, while its upper bound is O(L/k) in Lemma 6.4. However, the scale of difference in the absolute sense, even for the smallest k 20. Thus, while there may be gap between the empirical difference p(K+1) - p(( and the theoretical upper bound, the proposed method are expected to be stable practice.
21


Reviews:
[{"review_id": "QpxwzlL29G", "reviewer": "Reviewer_X4ks", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n- To the best of my knowledge, the approach is reasonably novel. The framework is general-purpose and the choice of Bezier simplexes seems well motivated. \n- The theoretical guarantees of the proposed method are a strong contribution. \n- Overall, the paper presents a reasonable method for a particular class of problems, and is very close to publication ready. The paper is quite clear and reasonably self-contained. \n\n# Weaknesses\n- The proposed approach is limited to simplicial problems. This is not a major weaknesses, as it is reasonable to restrict the set of problems of interest and exploit the structure of those problems. \n- The experimental evaluation seems (to a non-expert in this subfield) to be quite limited. I will leave it to the other reviewers to comment on whether the current set of experiments is comprehensive enough for publication. However, related to the last point: it would be interesting to experimentally characterize the performance on non-simplicial problems. Are such experiments possible?\n- As far as I can tell, there was no analysis of the scalability of the method. How many points are reasonable to define the hypersurface (in varying dimensions)? What is the computational complexity of the simplex fitting procedure? Are there any other factors impacting the scalability of the method over alternate approaches? \n- Could the gradient descent method be extended to multiple steps of gradient descent? If so, it could enable re-use of previous {x} in modeling, if this is useful. How would this impact the computation complexity of Bezier simplex parameter updates? Moreover, could this provide a bridge between sampling-based Pareto front computation and your approach?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents an approach for fitting Pareto frontier (parametric) surrogate functions in multi-objective optimization. In particular, it builds on Bezier simplexes with a simple algorithm that alternates between sampling points from a currently-fit simplex, optimizing (or otherwise updating) these points, and fitting a simplex to the (updated) points. The authors further discuss a version of this approach in which gradient descent is used to optimize individual points, and present an update algorithm for simplex parameters. The authors present a PAC stability analysis of the methods, which appears correct. The authors also include numerical experiments. I am not an expert in multi-objective optimization and so I do not know if the scope and scale of these experiments is reasonable for the subfield.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "As I mentioned previously, I believe this paper is quite close to ready for publication. Please address the questions in my weaknesses section above and, if possible, include the experiments mentioned and if reasonable.", "broader_impact_concerns": "This is a narrowly technical paper that I do not believe requires a broader impact statement.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "dTafXrry7P", "reviewer": "Reviewer_tznU", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n+ A new algorithm is built to transform a single-objective optimization \ninto a multi-objective optimization problem. \n\n+ Theoretical properties such as PAC stability of the algorithm is studied. \n\nWeaknesses: \n- Only a simple numerical simulation is performed to compare the proposed method. \nNo real world data is used for comparison. \n\n- It would be desirable if the stability property can be verified with the simulation or \nempirical results. \n\n- The results in Table 3 appear to show that the proposed method and the baseline have comparable accuracy \nand convergence speed. The advantage of the proposed method does not appear to be \nso much.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Overview: \n\nThis paper constructs an algorithm to convert a single objective optimization to a multiple \nobjective optimization and examine the PAC stability (and thus upper bound for the generalization gap).", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please address the weaknesses listed above. \n* verification of the stability empirically\n* experimental results on real world data\n* show clear advantages over baseline methods", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "56iAZp0Qft", "reviewer": "Reviewer_Df3Q", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n- Novel algorithm that provides significant conceptual and practical improvments\n- Novel generalization bounds for multi objective optimization algoorithms\n\n# Weaknesses\n## PAC stability and generalization\nThough the motivation for the PAC stability framework is to obtain generalization bounds, the effective obtention of such bounds for Alg. 2 is conditional on Alg. 2 being \"decomposable\", which does not seem trivial to me (cf remark after Thm. 6.5).\n\n## Numerical experiments\nThere is no mention of the computational overhead of the proposed method. Instead of performing a Bezier fit once like existing methods, this algorithm performs one at every iteration. It would be fairer to compare methods for a given computational budget, or at least report CPU times.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "- This work proposes a multi-objective optimization algorithm: it is the first method that directly operates on a parametrization of the Pareto set instead of working with a finite set of samples.\n- It then extends the stability of optimization algorithms to this setting and establishes generalization bounds for their algorithm.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- Discuss the \"decomposability\" issue.\n- Discuss the computational overhead.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This manuscript presents a novel approach to extend a single-objective optimization algorithm to the multi-objective setting by appealing to a tool known as the Bézier simplex model and provide an analysis of this approach in the PAC framework.

The paper received a positive reception during the initial review phase, with the reviewers only raising some minor questions and comments. The discussion phase served to clear up these minor issues, which were satisfactorily addressed by the authors in their comments and revision.

After the discussion / revision phase, the reviewers are unanimous in their support of the paper being accepted at TMLR. In particular, there is broad agreement that the claims made in the paper are supported by clear and convincing evidence and that the content of the paper would be of interest to a segment of TMLR's target audience.
Decision: Accept (poster)