Neural Evolution Strategy for Black-box Pareto Set Learning
Authors: Chengyu LU
Zhenhua Li
Xi Lin
Ji Cheng
Qingfu Zhang
Venue: NeurIPS 2025 poster ()
Keywords: evolution strategy
multi-objective optimization
black-box optimization
pareto set learning
Primary area: optimization

Abstract: Multi-objective optimization problems (MOPs) are prevalent in numerous real-world applications. Recently, Pareto Set Learning (PSL) has emerged as a powerful paradigm for solving MOPs. PSL can produce a neural network for modeling the set of all Pareto optimal solutions. However, applying PSL to black-box objectives, particularly those exhibiting non-separability, high dimensionality, and/or other complex properties, remains very challenging. To address this issue, we propose leveraging evolution strategies (ESs), a class of specialized black-box optimization algorithms, within the PSL paradigm. Traditional ESs capture the complex dimensional dependencies less efficiently, which can significantly hinder their performance in PSL. To tackle this issue, we suggest encapsulating the dependencies within a neural network, which is then trained using a novel gradient estimation method. The proposed method, termed Neural-ES, is evaluated using a bespoke benchmark suite for black-box PSL. Experimental comparisons with other methods demonstrate the efficiency of Neural-ES, underscoring its ability to learn the Pareto sets of challenging black-box MOPs.

Neural Evolution Strategy for Black-box Pareto Set Learning
Chengyu Lu 12, Zhenhua Lis.,,, Xi Ji Cheng!.?, Qingfu Zhang 1.2.
City University of Hong Kong, City University of Hong Kong Shenzhen Research Institute Nanjing University of Aeronautics and Astronautics
MIIT Key Laboratory of Pattern Analysis and Machine Intelligence
{chengyulu3-c, xi.lin, J. Cheng)@my cityu. edu. hk, zhenhua 1 Li@nuaa. .dd. cn qingfu. shangeeityy..dd.. edu.hk
Abstract
Multi-objective optimization problems (MOPs) are prevalent in numerous real- world applications. Recently, Pareto Set Learning (PSL) has emerged as a powerful paradigm for solving MOPs. PSL. can produce a neural network for modeling the set of all Pareto optimal solutions. However, applying PSL to black-box objectives, particularly those exhibiting non-separabilltyy high dimensionality, and/or other complex properties, remains very challenging. To address this issue, we propose leveraging evolution strategies (ESs), a class of specialized black- box optimization algorithms, within the PSL paradigm. Traditional ESs capture the complex dimensional dependencies less efficiently, which can significantly hinder their performance in PSL. To tackle this issue, we suggest encapsulating the dependencies within neural network, which is then trained using a novel gradient estimation method. The proposed method, termed Neural-ES, is evaluated using bespoke benchmark suite for black-box PSI.. Experimental comparisons with other methods demonstrate the efficiency of Neural-ES, underscoring its ability learn the Pareto sets of challenging black-box MOPs.
Introduction
Many rial-life applications need to consider multiple conflicting objectives simultaneously 133
Since these objectives can rarely be optimized by single solution simultaneously, commonly used approach is to find the Pareto set (PS) [4] which contains all solutions with different optimal trade- offs among the objectives. Traditional methods (5-)) approximate the PS using only a finite number of solutions, risking not aligning with users' trade-off preferences determined after optimization.
Recently, a novel paradigm called Pareto Set Learning (PSL) has attracted increasing attention [9] [5]. PSL formulates multi-objectivy optimization problem (MOP) 4 as a set learning task by decomposing the problem into an infinite number of single-objective subproblems and resolving the optimal solutions for all of them. In this process, a set model, typically neural network, learns the mapping from an input preference to its corresponding Pareto optimal solution. Consequently, the Pareto set is recovered in both the objective and solution spaces. Unlike traditional multi -bjei-ooeeectve methods, which only accommodate a limited number of predefined user preferences [16] 5 a trained PSL model can provide a tailored optimal solution for any valid decision-maker preference trade-off without requiring re-optimization from scratch.
Despite the rapid development of PSL. it remains largely unclear how to apply this paradigm to black-box MOP where both the analytic expression and derivatives the objectives are unknown. On
"Correspondence to: Qingfu Zhang and Zhenhua Li
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

one hand, black- box objectives are prevalent in practice 17 18 71, making it imperative to develop novel approaches for effectively tackling black-box MOPs. On the other hand, when exhibiting common properties such as non-separability and high dimensionality, black-hox objectíves can become highly challenging, rendering their resolution an open question. Moreover, since PSL. involves solving numerous subproblems that may themselves be black boxes, jointly optimizing them poses a significantly greater challenge than solving a single black-box optimization problem. Lin et al. [9] proposed integrating Gaussian process model into PSL for solving expensive MOPs. However, Gaussian process models scale poorly with problem dimensions and increased evaluation budgets, which may be less suited for more general black-box MOPs.
Evolution Strategy (ES) has remained one of the most successful black-box optimization algorithms for decades [9 21] An ES algorithm utilizes distribution model, typically multivariate Gaussian distribution, to search in complex fitness landscapes. Key elements that govern a competitive ES algorithm include an adaptive covariance matrix, the resemblance to gradient search, and sophisticated adaptation strategies. Compared to other zeroth-order optimization methods that sample solutions from isotropic distributions [22 23 ESs excel at dealing with variables' complex dimensional dependencies. a prevalent feature as well as a significant challenge in real-world optimization prob- lems 26 It is natural to integrate PSL with ES for black-box multi-objective optimization. On one hand, PSL decomposes an MOP into single-objective subproblems. aligning perfectly with the strengths of ES. On the other hand, ES's ' search efficiency paves the way for dealing with challenging black-box box objectives. However, traditional ES algorithms experience high computational complex- ity [27 28 and rely heavily adaptation strategies [1] These limitations particularly crucial and evident within PSI., as they can lead to failures in decoupling dimensional dependencies and significantly hinder the performance of ES.
26
To enhance the PSL paradigm for solving black-box MOPs, we propose a novel method called Neural- ES, which acts as the back-end of the set model. To effectively handle dimensional dependencíes, we propose rethinking the way ES algorithms decouple such dependencies. Rather than leveraging a covariance matrix, Neural-ES encapsulates the dependencies within neural network. To be more specific. we partition the high- digenniiensional design variable of a black-box objective into lower- dimensional segments. Each segment follows distinct multivariate Gaussian distribution. The network, which also acts as the set model, parameterizes the distribution conditioned on all previously sampled segments. To train the network, we estimate the gradient in novel way closely related to Neural-ES is well-suited for black-box tasks involving neural networks, contributions in threefold:
the natural gradient (32 33
with PSL serving as typical example. We summarize
We propose Neural-ES, novel algorithm that utilizes neural network to capture complex dimensional dependencies effectively, which enhances the versatility of the ES repertoire.
We establish a novel PSL formulation for black-box MOPs, which enables the use of efficient serivative-free optimization algorithms to address the problem efficiently.
We conduct extensive experiments using novel benchmark suite designed for evaluating black-box PSL problems. Results demonstrate that Neural-Er effectively addresses non- separable MOPs. To the best of our knowledge, are the first to scale the PSL paradigm to over 1,000 dimensions in black-box box scenarios.
2 Preliminaries & related works
2.1 Multi-ohjective optimization
continuous MOP (MOP) 4 is defined as:
minimize F(w)= = (.),
#EX
where F RM consists of M objectives X for all i 1,..,, M. many cases, the objectives conflict with each other, and single solution cannot optimizes all of them simultaneously. An alternative is to pursue the Pareto set, which consists of optimal solutions under all possible trade-offs among the objectives. Collectively, these solutions exclude any opportunity to improve one objective without compromising at least one other. We provide formal definitions of the fundamentals of MOPs in the following:
2

---

Definition 1 (Pareto dominance). Let x, x E X' be two solutions to (1. is said to dominate æ, denoted by æ x. f( S (z), vi = L, M, , and fz < f; (x) Ej = 1,..., A.....M.
* æ, Vo E X. Definition (Pareto optimality). A solution æ X is Pareto optimal
Definition (Pareto set & front). The of all the Pareto optimal solutions is called the Pareto set, denoted by PS {I E X É x / the image of PS in the objective space called the Pareto front, denoted by PF ((x) a E PS).
Multi-objective evolutionary algorithms (MOEAs) MOEAs have developed as one of the key methodologies for solving MOPs [66 5 6). Compared to other methodologies [34] 71. MOEAs
approximate the PS by population of diverse non-dominated solutions, without collapsing into single point. Several highly efficient MOEA variants have been developed [6 35 36). How- ever, MOEAs encounter longstanding issues, including inefficient Pareto ranking, misalignment with decision-makers preferences, poor scalability in relation to population size, and ineffective information sharing among subproblems.
2.2 Pareto set learning
PSL formulates MOP follows. First, given a preference vector (a),..., aM)) AM where AM = fo E | | a = J is the (M1--1ddmenioon simplex, and an aggregation function IR + IR conditioned on a, can be aggregated into a single-objective subproblem. Under some mild conditions, solving every possible subproblem conditioned on all possible a recovers the entire Pareto set. Then, an MOP in ①I can be formulated as:
minimize Ea-u((A))) [g oz)], where a = p)(a), BE
(2)
where θ is the parameter space, ((A denotes probability distribution (in this work, is set as uniform distribution) over AM. and Po AM - X is called a set model parameterized by 0. In the ideal case, (a ) = argmin g(æ; ox), Va E Ar, which means yo maps each preference to its corresponding Pareto optimal solution that minimizes a();
Recent developments Since its inception, PSL 9 has been rapidly applied to various problem types, including expensive [13 combinatorial (37], and smooth MOPs (38]. However, these methods face challenges with general black- box MOPs, particularly when the objectives are high- dimensional and exhibit complex dimensional dependencies. Lin et al. [15] contributed primary work called evolutionary-PSI (EPSL), where finite difference approach employed to train the set model. However, due the absence of comprehensive distribution model, most existing derivative-free PSL algorithms have been validated only on low-dimensional toy problems, leaving their performance challenging problems uncertain.
2.3 Evolution strategy
ES a prominent class of derivative-free optimization algorithms. Popular algorithms include natural evolution strategy (NES) [28] 20), the state-of-the-art covariance matrix adaptation (CMA-
ES) 19 27), and the also competitive but more time-efficient matrix adaptation evolution strategy
(MA-ES) 39 An ES algorithm first formulates general black-box optimization problem in question as:
minimize ((x) maximize -En(s [f())], øEX
where X' C Rd is the search space, 7 (xx; sp) is the probability density function, p.d.f., of smooth distribution x ~ π (sp) parameterized by and Φ is the parameter space of At each E X,, only its ff function value f(x) ) is available. Most commonly, π (5p) Nmm, C)) is multivariate Gaussian distribution where consists of the distribution mean m E R² and the covariance matrix C St (where i denotes the space of n-dimensional symmetric and positive definite matrices). In brief,
an ES optimizes of by 1) sampling random samples from IT (T()) and 2) effectively updating sp
Covariance matrix and dimensional dependencies widely acknowledged that the covariance matrix plays vital role in the performance of ES algorithm. Hansen et al. [24] systematically
3

---

showed that CMA-ES, aided by an adaptive covariance matrix, effectively decouples the interacting problem dimensions. Beyer [26] showed that the covariance matrix decouples the dimensional dependencies by effectively converging to the inverse of convex objective's Hessian. Both Aki- moto al. [40] and Glasmarchers el al. (28] demonstrated that the rank-j of CMA-ES amounts to the natural gradient ascent method. The ability to decouple dimensional dependencies positions NES and CMA-ES as second-order quasi-Newton methods rather than merely being zeroth-order finite-difference methods.
Adaptation strategies In addition to natural gradient, refined adaptation strategies also contribute to the competitiveness of ESs, as exemplified by cumulative step-length control and evolution paths [19 41 An evolution path accumulates the movement trajectories the distribution over time, thereby exploiting promising search direction. Evidence suggests that ESs leveraging both natural gradient and adaptation strategies significantly outperform those that rely solely on natural
gradient [28 29]
3 Black-box Pareto set learning
Preference
Set model (
Distribution a(:
Solution space X
Objective
a
F(O
a
f
Ay
Sampling x ~ m(:
5t m) D.+
ad A 2xp 88
Update 0
Figure I: Black-box PSL formulation. The blue arrows indicate the forward process, and the yellow arrow indicates the gradient pack-propagation: the shadowed area is where an ES participates and is referred as the buck-end of the model.
3.1 Problem formulation
black-box MOP refers to those in that provide no information beyond the objective values of solution. In this case, we formulate through the reparameterizatioo trick:
maximize [(0) Ea-U(s Ba(erp)] (=)(;; g(x; ox)]. where = (a). BE
The above formulation 4 distinguishes itself from in the nested expectation. Instead outputting directly, the set model (e ((x) first maps preference vector to distribution π ((p)) parameterized by +p (a)). and then samples solutions from the distribution. To train the set model, an ES is employed to sample solutions and estimate gradients. When 0 is perfectly learned, regarding an arbitrary a, a random æ π (x))(ax)) converges in probability to the minimum of g(:
To further illustrate this process, we outline the paradigm of black-hox PSI. in Fig. First, a random preference vector a U(AM) is sampled and inputted into the set model (p.. The model subsequently outputs a parameter vector E ④, which parameterizes distribution ग ((;; sp). A random solution æ ((p) then sampled from this distribution and evaluated the specific
g(::a). a). Here, we use the PBI aggregation function phi [8] as an example. When we have sampled enough as, and for each there have been several samples æs, can estimate the gradient and back- pook-propagate it to 0 through the chain rule.
Key challenges
The formulation introduces several challenges for ES. First, the high computational complexity of the covariance matrix significantly impedes an ES's performance. It takes ((n²) space complexity of the set model to generate a covariance matrix, and ( time complexity or even more to compute the gradient, which will soon become prohibitive as n increases. Even worse, since each distribution yp corresponds to different can be seen from Fig. the large number of decomposed subproblems

---

in PSL necessitates an equally large number of covariance matrices, which drastically increases the learning burden of the model.
Second, it is difficult, not impossible, to apply sophisticated adaptation strategies in PSL. problems. These strategies are known accelerate convergence by counteracting the slow-changing rate suppressed by the Kullback-Leibler -KLe-diveeeence [26]. Despite some recent findings [42]. the connection between these strategies and gradient search remains unclear, which limits their use in training the set model. Furthermore, as there are numerous decomposed subproblems in a PSL problem, infeasible to track a different evolution path for each subproblem.
The above characteristics black-box PSL problem make the direct integration of ESs rarely yield competitive performance, posing new challenges for ESs to learn in typical adaptation-fre -feee scenario.
4
Neural evolution strategy
4.1 Basic idea
Consider a general single-objective problem: maximizeg [o( = (m(pp) (f where sp < a is problem specific parameter, and æ E RR. Assuming n = p x d where and d are positive integers, z can be partitioned into p segments each of length d, as follows:
x X1 Xdy [(I (il-))--::.- Xids
Xn d+1 Xp
(5)
By conditional probability, TT (p) a (æ; uoo (a)) can be decomposed into the product of chain of p.d.fs each concerning a different x;, and conditioned on all the previous To, E. .... I;-1:
x(;; (p)))) II No a zo, 2 04-1
N (x); 4e (2):--))))
(6)
ma.
With a slight abuse of symbols, we redefine Po U._ IR' - R2 x $4 hereafter, which means (e maps variable-length inputs to the pa- rameters of d-dimensional Gaussian. More specifically, input Ø0.-1 consists of all the segments Vo, Ti..., x;-1 preceding I;, and the output spe (0);--1) consists of the mean mo(#o.i-)) E Rd and the covariance matrix Ce( E se_.. Zo is just an alias of a. for notation consistency. (6) directly models
Xa
the dependency of æ; on #0:i-1 as the mapping so. practice, sp is neural network. Fig. 2
Figure 2: An illustration of using (6) to sample random solution T(pp((a)).
illustrates the idea.
As the Gaussian is smooth w.r.f. sp, the network
can be trained by gradient ascent. A conceptual algorithm is inspired, which iteratively 1) sample (: (p) independently, and 2) estimate the gradient using these samples and back-propagate to 0 until convergence.
solutions (a> } from (
4.2 Implementation
4.2.1
Rank-one covariance matrix
For any E [p]. full covariance matrix Co (mo:-11 1) of O(d²) degrees of freedom necessitates (dd) output neurons to generate it and ((d)) time complexity to sample solutions as well as to estimate its gradient. To reduce the overall complexity to (nn), the above must be reduced to ((d).
As is common trick ES variants for dighedimensional tasks 43 covariance matrix by rank-one model. More specifically, we define a step-size do (Toa-1) > 0 and a component to (1)x-1) 1) E Rdd, then Co (x can be parameterized by:
29
30),
we can replace a full
Co(xo--1) (x) 1) (I vo (0 1) (xo
(7)
5

---

Of the above, de controls the isotropic variance of the samples, and is expected to identify promising search direction. Although the number of free parameters in Cp decreases, I (Tp)) suffices to decouple complex dimensional dependencies when is large enough. By this means, the space complexity of 0's output layer reduced (d).
Following 7 solution is practically generated by:
2i mo 10x--1 + 0o 00( + rvg (2oo--11)
(8)
where 2 ~ N0 I), E R., and (0, 1). After iterations, an ~ no(: a)) is recovered by stacking the 31, Tp. The time complexity of (8) is (d))
4.2.2 Natural gradient
The gradient of 0 can be derived as
θ (OInN (Poxxxx))))
(O)
Σ Ba(z:p)) i=1
f(x)
00
d (To:i-1)
However, the partial derivative IInNV(;;;; (m (0o: rarely yields stable conver- gence behaviors of the mean mo 1)) because the shrinking of Co (1) in turns stretches the gradient of mo(Z0o:-1), preventing it from staying in stationary points. See (20 26 as well as
Section [B2 for more details.
To regulate the gradient, a remedy is to multiply each partial derivative in ( by the inverse of the Fisher information matrix (a.k.a. the "Fisher") of yp: (00. which yields:
à [(x)
θ In (p (
)]
(J.
Σ E (axp))
(((p.))
(10)
a0
aupa (205-1)
where F(Pp:) = E N(ta/se.) VV. In, pP: In P.).. The interpretation is that the Fisher penalizes drastic changes of the search distribution, counteracts the effect of the shrinking covariance, and thus guarantees stable convergence. More intrinsically, the following connects (10) to the natural gradient, which incorporates second-orde information and is the steepest ascent direction on Riemannian manifolds [33]
Lemma Consider distribution ((p) with fixed parameters &p = (4T .........9 ,T
".. and ran- (p.) / in [0 is an dom sample a ~ T(p). The partial derivative Σ ["=,)((())) 10 Wn
unbiased estimator for the natural gradient of the log-tikelihood In (αx; sp) regarding 4p
have shown that unscaled objective values often lead to degenerated search distributions and premature convergence, due to the irregularity of ʃ. ES algorithms alleviate this issue by replacing (ax) with rank-dependent weight values. Specifically, consider set of objective values {f(x)) }} arranged descending order, which means f(x') 2 2 (x**, we replace each f(x) by wy where wh 2 ... 2 w).- In this paper, we define my = In (1 + Inj.
Furthermore, empirical and theoretical evidence 20 25
26
According Akimoto el al. [4)] the Fisher F(4::) can be absorbed into the gradient without explicitly computing its inverse. Eventually, the gradient estimator based on a set of solutions { J==1 can be derived as:
J(miy)) ((m m.;)+
Rol 112 + - od ||: [1 Ooiz
(0) ~J
00 , (11)
o o (w )" (10 )y + || - 20 (1 + | 2))),,))))))))))))))
(m where maj- = mg (ar.i-11) 1), oij = o |
Vig = u (o( Vij = m:j. J(-) computes
the d x (0) dimensional Jacobian matrix of an input vector w.r.t. 0, and || computes the norm.
We leave the tedious derivation steps in Section B.4 Each term irrelevant to 0 (i.e., those except
J (m.;), /0, and J ((:;))) yields (dd time complexity.
4.3 The algorithm

---

Algorithm 1: Neural-ES for black-box PSL
Input :A black-box MOP F(x)
We present Neural-ES in Alg. During each epoch, i.e, lines 2 8, the algorithm optimizes N randomly sampled subproblems g(;aa), E [N]. For each subproblem, lines 4 to 6 sample À random solutions from ((po(a)))) and line 7 computes the gradient regarding this subproblem. Lines 4 to 7 cor-4 respond to an iteration of typical ES algo-5 rithm. Line 8 estimates the gradient w.E.t. the PSL problem in Eo- VJ.+ (0)
Output: :A set model ye
while The termination is not triggered do for k 1 to N do
Sample preference at ~ u((m),
set 22g = a=: for j = 1 do
for i to p po Construct a. based
4 which reads ( E (O).
Compute V.a (0) based on
where f(-) == =²(;a²);
Eventually, (O) is back-propagated to 0, where 7) denotes the learning rate.
0 0 +
N Σ k=1 VJ.a (0);
5
Experiments
5.1
Configurations
Black-Box Pareto Set Learning (BBO-PSL) suite To evaluate the performance of PSL algorithms, a test problem should meet a few fundamental requirements, including 1) closed-form PF and PS, 2) a complex shape of the PS manifold, 3) scalability arbitrarily high dimensions, and 4) non- separability with a wide spectrum of the Hessian. However, most existing problem suites 444 either implement only few of these requirements
include other interfering features that fall outside the scope of this paper. To this end, we present the novel BBO-PSL suite, which comprises eight instances (FI to F8) and satisfies all the aforementioned requirements. Among them, Fl to F4 are bi-objective instances, and F5 to F8 are tri-objective.
Compared algorithms We compare Neural-RS 2 to the state-of-the-art and canonical EPSI. [5]. and configure Neural-ES by d = [n/6], of the most effective methodologies for obtaining set of non-dominated solutions, we consider two state-of-the-art MOEAs specifically designed for high-dimensional MOPs, namely LMO-CSO (35] and IM-MOEA/D [36]. Among the many MOEAs, LMO-CSO is favored for its superior handling of dimensional dependencies in complex and high-dimensional landscapes, while IM-MOEA/D makes notable effort to elevate the canonical MOEA/D g to high dimensions by incorporating the popular inverse modeling technique [48). Compared to most others, they are more lightweight with lower computational overhead, indicating better scalability to larger population sizes. Moreover, we synthesize an algorithm called EPSL-RI, by setting p in Neural-ES,, which represents the simplest way to integrate PSI. with an ESS
[5n/6]], and À = 2 + In]. Additionally,
MOEAs widely recognized
Metrics The Inverted Generational Distance (IGD) 49 is used the performance metric, with a smaller value indicating better perfor- mance. Each algorithm is repeated 21 times on each BBO-PSL instance using different ran- dom seeds. For the three PSL algorithms, we evaluate their testing performance. To do so, we input large number of preferences into trained set model, which yields the same num- ber of distributions; from each distribution, we randomly sample a solution, the collection of which form the testing solution set; finally, we measure the IGD value of the set and refer it as the testing IGD value. Meanwhile, we track the convergence of the PSL algorithms during training, in terms of the IGD values acquired by the best solutions far. The MOEAs directly optimize the
EPSL
EPSL-RI
Neural-Ee ES
10
12
14
pettrtions (x10D)
Figure 3: Convergence of the PSI. algorithms on F1-128D.
*Code of Neural ES and the BBO-PSL suite is available in https::/gitte com/ coon/channlerO)) Neural-Evolution-SS5
7

---

same number of solutions, and the IGD value of the final population is measured. To examine the scalability of the algorithms with increasing dimensions, we vary n from 32 1024.
5.2 Results and discussions
10
Rgart
(a) Fl
(b) F2
(c) F3
(d) F4
kgam
(e) F5
(f) F6
(g) F7
(h) F8
Figure 4: Average testing IGD values w.r. problem dimensions, where the red curve with circles refers to Neural-ES (ours), the blue curve with diamonds refers EPSL-R1, the green curve with squares refers to EPSL, the orange curve with upside-down triangles refers to LMO-CSO, and finally the black curve with triangles refers to IM-MOEA/D. Their detailed values, along with the Mann-Whitner Whitney U test, can be found in Table 3
According to Fig. Neural-ES outperforms all the other algorithms in all eight instances by large margin. Although the MOEAs are specialized for high-dimensional MOPs, they struggle to scale efficiently when users require numerous trade-offs among objectives. This drawback arises from their substantial demand on evaluation budgets and limited generalization ability. In contrast, the PSL paradigm performs well during surges ones emerge after optimization begins. The PSL algorithms, and Neural-ES in particular, are more robust against increasing problem dimensions. All these attributed to the efficiency of evolution strategies and the generalization ability of the set model.
demand, especially when demands shift or
PS EPSI EPSI-AI Moural
Neural ES
EPSL EPSL R1 Neural- PP
EPSL EPSL PF
(a) Pareto sets of
(b) Pareto sets of F2
(c) Pareto fronts F1
(d) Pareto fronts of F2
Figure 5: Testing solution sets of the algorithms, F1-512D and F2-512D.
Compared to EPSL and EPSL-RI, Neural-ES reduces the IGD values by effectively decoupling the dimensional dependencies. Although EPSL-RI also handles the dependencies, it performs slightly better than EPSL in just few cases, yet remains similar overall. And its advantage over EPSL dissipates
increases. This verifies that simply plugging ES algorithm into adaptation- free scenarios, such as black-box PSL, does not yield a competitive algorithm. contrast, the novel treatment of dimensional dependencies positions Neural-ES as highly efficient alternative in adaptation- free scenarios.
In the remainder, we limit our focus to the PSL paradigm and utilize the limited space for in-depth discussions. addition to the testing performance, we plot in Fig. the training-time convergence of the three PSL algorithms. Neural-ES establishes performance lead over EPSL and EPSL-RI

---

P5 215L-11
25
Neural-BS
(a) Pareto set by EPSL (b) Pareto set by EPSL-RI (c) Pareto set by Neural-ES
(d) Pareto fronts
Figure 6: Testing solution sets of the algorithms, on F5-128D. For the ease of visualization, we only demonstrate the non-dominated solutions in ((d) Therefore, the results of EPSL-RI are unreported because they were dominated.
during the early stages, and subsequently extends the performance gap throughout the search process. The results suggest that Neural-ES efficiently learns the complex dimensional dependencies from the outset and exploits this advantage to significantly reduce the optimization difficulty.
The ability to approximate Pareto manifolds during testing is a unique feature of PSL as well as its primary priority. We plot in Figs. and some example approximation results. For the bi-objective Fl and F2, in Figs. 5(a) and 5(b) the approximate PSs by Neural-ES align more closely with the
ground-truth and spread much wider than EPSL and EPSL-RI. It turns out that the approximate PFs
by Neural-ES dominate those by the other two almost everywhere, as seen in Figs. 5(c) and 5(d)
Figs. ((a)) to ((c) depiet the approximate PSs in the first three dimensions of the search space (i.e., the J1 - 22 Ø3 subspace the x-space). these figures, the gray surface represents the ground-trun truth PS. while the colored points indicate the sampled solutions obtained by different algorithms during testing. In Fig. [[a]] solutions sampled by EPSL appear scattered and poorly aligned with preference vectors, because 1) the finite difference approach does not minimize the uncertainty of distribution, and 2) EPSL does not learn the dimensional dependencies. Fig. (6d)) depicts the ground-truth PF (i.e., surface with gray contours) and the algorithm approximations in the objective space (i.e., the F-space), which provides a more rigorous performance assessment because directly evaluates solution quality in the objective space and reveals dominance relationships that are not apparent in the æ-space. While Neural-ES and EPSL-RI appear similar in decision space (Figs.]6(()) to ((C)). Neural-ES clearly outperforms in objective space, with most EPSL-RI solutions being dominated by Neural-ES solutions.
5.3 Real-world applications
Trajectory planning We consider a tri-objective trajec- tory planning task modified from the popular rover prob- lem [50 34]. The goal is to determine a trajectory that minimizes 1) the cost during navigation, e.g., collisions with obstacles, 2) the distance from the endpoint to an intended target location, and 3) the total trajectory length, e.g., for saving energy. The objectives conflict with each other because, for instance, minimizing the distance may lead to increased frequency of collisions. We plan the tra- jectory by designing a total of 50 (x.y) coordinates, which means n = 100. Their approximate PFs are depicted in Fig. 7 Measured a reference point of (5.2, Neural- ES, EPSLL-R1, and EPSL achieve an average hypervolume of 8.3440..1155 5.122, and 4. 4.186 709. respec- tively, with higher values indicating better performance, Additionally, for each algorithm, we collect the optimal solutions for 300 evenly distributed subproblems. 299 out of 300 solutions from EPSL were dominated by those from Neural-ES, while 194 from EPSL-RI were also dom-
EPSL EPSL-R1 Neural-ES
diafance
fy cosk
Figure 7: Approximate PFs of the trajec- tory planning problem.

---

inated by those from Neural-ES. The results demonstrate the superiority of Neural-ES in balancing the cost, distance, and total length of the trajectories.
Multi-objective Unmanned Aerial Vehicle (MO-UAV) Navigation We
consider a set of real-world multi- objective MO-UAV navigation prob-
Table I: Average HV values and standard deviations, on the MO-UAV problems, with larger values the better.
lems, which are featured in the Meta-
Case
n
EPSL
EPSL-RI
Neural-ES
BBO-v2 library (5I This problem consists of three competing objectives: 1) the distance from the UAV's final location to the destination, penalized by the trajectory's ss non-smoothness, 2) the risk collision with threats, such as ground-to-air missiles, and 3) altitude penalty maintain safe and cost-efficient flight. By varying the control parameters of the terrain, we have developed three instances: easy, moderate, and hard, which require a UAV trajectory of 10, 20, and 40 waypoints, respectively. Because each point comprises H, y. and coordinates, the overall problem dimensions are 30, 60, and 120, respectively. Complex dimensional dependencies arise from the spatial and temporal correlations between waypoints, making navigation highly challenging. After 21 independent runs, we report the brief results in Table [ Neural-FS demonstrates its superiority in practice with the highest HV values across all difficulty levels.
Easy
30
2.47
1.58
3.40.0.25
Moderate
60
2.4040.37
1.52=0.27
3.44
Hard
120
0.72+0.20
0.84=0.14
2.44
6 Conclusion, limitation, and future work
Conclusions In this work, we have made an initial effort to solve novel yet important problem: black-box Pareto set learning (PSI.). The goal is to learn the entire Pareto set of black-box multi- objective function. We have developed novel formulation, which facilitates the use of powerful black-box optimization algorithms, particularly evolution strategies. To address the challenges posed by black-box PSL, we have developed novel mechanism to decouple the complex dimensional dependencies of the objectives. The proposed Neural-ES excels at learning challenging Pareto sets, successfully scaling the PSL paradigm to over 1000 dimensions with competitive performance.
Limitation and future work While this work focuses on unconstrained MOPs, constraints often significantly increase optimization difficulties, making constrained MOPs critical area to address in PSL. Developing novel and advanced methods to extend PSL to constrained MOPs is promising yet challenging direction for future research. Besides, Ye et al. recently suggest that when reparameterizing an aggregation function g ax as (=))))))) [g [ ). the optimization may only guarantee an approximate Pareto stationary point rather than exact one, in regardless of the specific distribution æ. It is very interesting to further investigate what led to such a mismatch. Nonetheless,
believe that the formulation in (4) remains of the best options available.
Acknowledgments
The work described in this paper was supported by the Research Grants Council of the Hong Kong Special Administrative Region. China [GRF Project No. CityU-11215723) and the Natural Science Foundation of China (Project No: 62276223). Zhenhua Li is supported in part by the Fundamental Research Funds for the Central Universities under NJ202433...
http::://gthhbbb.oo/o///bbbbbbtb
10

---

References
[1] K. M. Jablonka, G. M. Jothiappan, S. Wang. B. Smit, and B. Yoo, "Bias free multiobjective active learning for materials design and discovery." Nature communications, vol. 12, no. I. p. 2312, 2021.
O. Sener and Koltun, "Multi-task learning as multi-objective optimization." Advances in neural informa- tion processing systems, vol. 31, 2018.
[3] G. P. Rangaiah, Multi-objective optimization: techniques and applications in chemical engineering, vol. 5. world scientific, 2016.
[4] K. Miettinen, Nonlinear multiobjectivee optimization, vol. 12. Springer Science & Business Media, 1999.
[5] K. Li. "A survey of multi-objective evolutionary algorithm based on decomposition: Past and future," IEEE Transactions on Evolutionary Computation, 2024.
[6] T. Huang, S. Wang, and K. Li. "Direct preference-based evolutionary multi-objective optimization with dueling bandits," Advances in Neural Information Processing Systems, vol. 37, pp. 122206 172206-122258 122258, 2024.
[7] F. Ye, Y. Lyu, Wang, Y. Zhang, and Tsang, "Adaptive stochastic gradient algorithu for black-box oulti-objective learning," in The Twelfth International Conference on Learning Representations. 2024.
[8] Q. Zhang and H. Li, "MOEA/D: A multiobjectivv evolutionary algorithm based on decomposition." JEEE Transactions on Evolutionary Computation, vol. 11, no. 6, pp. 712-731, 2007.
[9] X. Lin, Z. Yang, X. Zhang, and Q Zhang, "Pareto set learning for expensive multi cuutti-obeective optimization." in Advances in neural information processing systems, vol. 35, pp. 19231-19247, 2022.
[10] E. Liu, Y.-C. Wu, X. Huang. C. Gao, R.-J. Wang, K. Xue, and C. Qian, "Pareto learning for multi- objective reinforcement lcarning," in Proceedings of the AAAI Conference on Artificial Intelligence, 2025.
[11] Y. Lu, B. Li, and A. Zhou, "Are you concerned about limited function evaluations: Data-augmented pareto set learning for expensive multi-objectto optimization." in Proceedings of the AAAI Conference on Artificial Inteliience, pp. 14202-14210, 2024.
[12] R. Yc, L. Chen, J. Zhang, and H. Ishibuchi, "Evolutionary preference sampling for parcto set learning," in Proceedings of the Genetic and Evolutionary Computation Conference, pp. 630-66 030-638,770 2024.
[13] M.-D. Nguyen, P. M. Dinh, Q.-H. Nguyen, L. Hoang, and D. D. Le, "Improving pareto learning for expensive multi-ohjective optimization via stein variational hypemetworks." in Proceedings of the AAAI Conference on Artificial Intelligence, pp. 19677-19685, 2025.
[14] Y. Liu, J. Yang, X. Ren, Z. Xinyi, Y. Liu, B. Song. X. Zeng, and H. Ishibuchi, "Multt-objecttive molccular design through learning latent pareto set," in Proceedings of the AAAI Conference on Artificial Intelligence, pp. 19006-19014, 2025.
[15] X. Lin, x. 7hang, 7. Yang, and Q. Zhang, "Dealing with structure constraints in evolutionary pareto set learning." IEEE Transactions on Evolutionary Computation, 2025.
[16] G. Lárraga and K. Miettinen, "Survey of interactive evolutionary decomposition- multiobjective optimization methods." Evolutionary Computation, pp. 1-39, 2025.
[17] Y. Shi, K. Xue, S. Lei, and C. Qian, "Macro placement by wire-mask-guided black-box optimization," Advances in Neural Information Processing Systems, vol. 36, pp. 6825-68668 6843,, 2023.
[18] P. N. Williams and K. Li, "Black-box sparse adversarial attack via multi-objective optimisation." in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 12291-12301, 2023.
[19] N. Hansen and A. Ostermeier, "Completely derandomized self self-adaptation in evolution strategies," Evol. Comput., vol. 9, no. pp. 159-195. 2001.
[20] D. Wierstra, Schaul. T. Glasmachers, Y. Sun, Peters, J. Schmidhuber. "Natural evolution strategies." The Journal of Machine Learning Research, vol. 15, no. 1, pp. 949-980, 2014.
[21] R. T. Lange, T. Schaul, Y. Chen, T. Zahavy, V. Dalibard, C. Lu, S. Singh, and S. Flennerhag, "Discovering evolution strategies via meta- meta-black-Ba black optimization," in The Eleventh International Conference on Learning Representations, 2023.
11

---

[22] O. Li, J. Harrison, J. Sohl-Dickstein, V. Smith, and L. Metz. "Variance-reduced gradient estimation via noise-reuse in online evolution strategies." Advances in Neural Information Processing Systems, vol. 36, pp. 45489-45501, 2023.
[23] P. Vicol, L. Metz, and J. Sohl-Dickstein "Unbiased gradient estimation unrolled computation graphs with persistent evolution strategies," in International Conference on Machine Learning, pp. 10553-10563, PMI.R, 2021.
[24] N. Hansen, R. Ros, N. Mauny, M. Schoenauer, and A. Auger. "Impacts of invariance in search: When CMA-ES and PSO face ill-conditioned and non-separable problems." Applied Soft Computing, vol. 11, no. 8, pp. 3755-5769, 2011.
[25] Y. Ollivier, I., Amold, A. Auger, and N. Hansen, - optimization algorithms: A unifying picture via invariance principles." Journal of Machine Learning Research, vol. 18, no. 18, pp. 1-65, 2017.
[26] H.-G. Beyer, "Convergence analysis of evolutionary algorithms that are based on the paradigm of informa- tion geometry." Evolutionary Computation, vol. 22, no. 4. pp. 079-709, 2014.
[27] N. Hansen, S. D. Müller, and P. Koumoutsakos, "Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)," Evolutionary computation, vol. 11, no. 1, pp. 1-18, 2003.
[28] T. Glasmachers, T. Schaul, S. Yi, D. Wierstra, and J. Schmidhuber, "Exponential natural evolution strategies." in Proceedings of the 12th annual conference on Genetic and evolutionary computation, pp. 393 400, 2010.
[29] Y. Akimoto, Auger, and N. Hansen, "Comparison- natural gradient optimization in high dimension," " in Proceedings of the 2014 Annual Conference on Genetic and Evolutionary Computation, pp. 373-380, 2014.
[30] 7. Li and Q. Zhang, "A simple yet efficient evolution strategy for large-scale black-box optimization," "
IEEE Transactions on Evolutionary Computation, vol. 22. no. 5. pp. 637-646, 2017.
[31] I. Loshchilov, T. Glasmachers, and H.-G. Beyer, "Large scale black-box optimization by limited-memory matrix adaptation," IEEE Transactions on Evolutionary Computation, vol. 23, no. 2, pp. 353 358, 2018. [33] J. Martens, "New insights and perspectives the natural gradient method," Journal of Machine Learning Research, vol. 21, no. 146, pp. 1-76, 2020.
[32] S.-L. Amari, "Natural gradient works efficiently in learning," Neural computation, vol. 10, no. 2, pp. 251- 276, 1998.
[34] S. Daulton, D. Friksson, M. Balandat, and E. Bakshy, "Multi-ohjective bayesian optimization over high- dimensional search spaces," in Uncertainty in Artificial Intelligence, pp. 507-517, PMLR, 2022.
[35] Y. Tian, X. Zheng, X. Zhang, and Y. Jin, "Efficient large-scale multiobjective optimization based on a competitive swarm optimizer," IEEE Transactions on Cybernetics, vol. 50, no. 8, pp. 3696-3708, 2020.
[36] L. R. Farias and A. F. Araújo, "IM-MOEA/D: An inverse modeling multi-objective evolutionary algorithm based on decomposition." in 2021 IEEE international conference on systems. man, and cybernetics (SMC). pp. 6664467, IEEE, 2021.
[37] X. Lin, Z. Yang, and Q. Zhang, "Pareto set learning for neural multi-objective combinatorial optimization," in International Conference on Learning Representation, 2022.
[38] X. Zhang, X. Lin, B. Xue, Y. Chen, and Q. Zhang, "Hypervolume maximization: A geometric view of pareto set learning"" " Advances in Neural Information Processing Systems, vol. 36, pp. 38902-38929, 2023.
[39] H.-G. Beyer and B. Sendhoff, "Simplify your covariance matrix adaptation evolution strategy," IEEE Transactions on Evolutionary Computation, vol. 21, no. 5, pp. 746-759, 2017.
[40] Y. Akimoto, Y. Nagata, I. Ono, and S. Kobayashi, "Bidirectional relation between ema evolution strategies and natural evolution strategies." in Parallel Problem Solving.roum Nature, PPSN XI: 11th International Conference, pp. 153-163, Springer, 2010.
[41] N. Hansen and A. Auger, "Principled design of continuous stochastic search: From theory to practice," in Theory and principled methods for the design of metaheuristics, pp. 145-180, Springer, 2013.
12

---

[42] R. Harmano, S. Shirakawa, and M. Nomura, "Natural gradient interpretation of rank-one update in CMA-ES," in International Conference on Parallel Problem Solving from Nature, pp. 252-267, Springer, 2024.
[43] Y. Sun, T. Schaul, F. Gomez, and J. Schmidhuber, "A linear time natural evolution strategy for non- separable functions," in Proceedings the 15th anmual conference companion on Genetic and evolutionary computation, pp. 61 61-62, 2013.
[44] Q. Zhang, A. Zhou. Zhao, P. N. Suganthan, W. Liu, Tiwari, et al., "Multiobjective optimization test instances for the cec 2009 special session and competition," tech. rep.. Colchester, UK, 2008.
[45] D. Brockhoff, A. Auger, N. Hansen, and T. Tušar, "Using well-understode single-objective functions in multiobjective black-box optimization test suites," Evolutionary Computation, vol. 30, no. 2. pp. 165-193, 2022.
[46] R. Cheng, Y. Jin, M. Olhofer, et al., "Test problems for large-scale multiohjective and many-objective optimization." IEEE transactions cybermetics, vol. 47, 12, pp. 4108 4121, 2016.
[47] S. Zapotecas- Martinez, C. A. C. Cocllo, H. E. Aguirre, and K. Tanaka, "Challenging test problems for multi-and obje-tbeective optimization." Swarm and Evolutionary Computation, vol. 81, 101350, 2023.
[48] R. Cheng, Y. Jin, K. Narukawa, and B. Sendhoff, "A multiobjectivve evolutionary algorithm using gaussian process-based inverse modeling," " IEEE Transactions on Evolutionary Computation, vol. 19, no. 6, pp. 838- 856, 2015.
[49] J. G. Falcón-Cardona and C. A. C. Cocllo, "Indicator-based multi-objective evolutionary algorithms: A comprehensive survey"" ACM Computing Surveys (CSUR), vol. 53, no. pp. 1-35, 2020.
[50] Z. Wang, C. Gehring, P. Kohli, and S. Jegelka, "Batched large-scale bayesian optimization in high- dimensional spaces," in International Conference on Artificial Intelligence and Statistics, pp. 745-754, PMLR, 2018.
[51] Z. Ma, Y.-J. Gong, H. Guo, W. Qiu, S. Ma, H. Lian, J. Zhan, K. Chen, C. Wang, Z. Huang, et al., "Metabox- v2: A unified benchmark platform for meta-black-bon optimization," arXiv preprint arXiv:2505.17745 2025.
[52] K. Gao and O. Sener, "Generalizing gaussian smoothing for random search," International Conference on Machine Learning, pp. 7077-7101, PMLR, 2022.
[53] S. Boyd and L. Vandenberghe, Convex optimization. Cambridge university press, 2004.
[54] S. Yi, D. Wierstra, T. Schaul, and J. Schmidhuber, "Stochastic search using the natural gradient," in Proceedings of the 26th annual international conference on machine learning, pp. 1161-1168, 2009.
[55] Z. Li and Q. Zhang, "What does the evolution path learn in ema-es?," in International Conference on Parallel Problem Solving from Nature, pp. 751-760, Springer, 2016.
[56] Y. Akimoto and N. Hansen, "Diagonal accelcration for covariance matrix adaptation evolution stratogics," Evolationary computation, vol. 28, no. 3. pp. 435, 2020.
[57] Y. Akimoto and N. Hansen, "Projection- restricted covariance matrix adaptation for high dimension." in Proceedings of the Genetic and Evolutionary Computation Conference 2016, pp. 197-204, 2016.
[58] Y. Lyu and I. W. Tsang, "Black-box optimizer with stochastic implicit natural gradient," in Joins European conference on machine learning and knowledge discovery in datobases, pp. 217-232, Springer, 2021.
[59] Y. Lyu, "Fast rank-1 lattice targeted sampling for black-box optimization," Advances in Neural Information Proce Processing Systems, vol. 36, pp. 10494 10515, 2023.
[60] B. L., J. K. Tang, and X. Yao. "Manyvobeeetive evolutionary algorithms: A survey," ACM Computing Surveys (CSUR), vol. 48, no. I, pp. 1-35, 2015.
[61] D. A. Harville, "Matrix algebra from a statistician's perspective." 1998.
[62] K. Deb. L. Thiele, Laumanns, and E. Zitzler, "Scalable test problems for evolutionary multiobjective optimization," " in Evolutionary multiobjective optimization: theoretical advances and applications, pp. 105- 145, Springer, 2005.
[63] H. L. and Q. Zhang, "Multiobjective optimization problems with complicated pareto sets, moca/d and nsga-ii," IEEE transactions on evolutionary computation, vol. 13, no. 2. pp. 284-302, 2008
[64] Y. Tian, R. Cheng. X. Zhang, and Y. Jin, "PlaEMM: A MATLAB platform for evolutionary objei-obeective optimization," IEEE Computational Intelligence Magazine, vol. 12, no. 4, pp. 73 87, 2017.
13

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction clearly state the claims, including the contributions made and limitations. We carefully exclude any over-claim our contributions.
Guidelines:
The answer NA means that abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match thcoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
This paper discusses constrained multiobjective optimization problems, a prevalent area that remains unsolved.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate" "Limitations" section in their paper.
paper should point out any strong assumptions and how robust results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model well-specification,
asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested a few datasets with few runs. In general, empirical results often depend implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to text system might not be used reliably to provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
appplicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerming limitations.
3. Theory assumptions and proofs
For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
Question:
Answer: [Yes]
We provide detailed proofs and derivations in the supplementary materials. All theorems, formulas, and proofs in the paper numbered and cross-referenced. All assumptions are clearly stated or referenced in the statement of any theorems.
Guidelines:
The answer NA means that the paper does not include theoretical results.
14

---

All the theorems, formulas, and proofs in the paper should be numbered and cross cross-reference..
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can cither appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by
formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: The steps of our proposed algorithms are deseribed step by step. The code will be open-sourced upon article acceptance for reproducibility experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, contribution novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data often good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (c.g., in the case of a large language model), relcasing model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility. which may depend on the nature of the contribution. For example
(a) If the contribution is primarily new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g.. a large language model), then there should either be a way to access this model for reproducing the results a way to reproduce the model (e.g.. with open- dataset instructions for how construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility. In the case of closed source models, it may be that access to the model is limited in some way (e.g.. registered users), hut it should be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: The code will be open-sourced upon article acceptance.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines for more details.
http:////ipps.ccppppbeee
guides/Cood Sataesss/oooSSuuisiiiiioooooo
While we encourage the release of code and data, we understand that this might not be possible,
so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
15

---

The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines htttps://nipa.cc/ppun for more details.
Buidee//DodeSuumiissioooPPo
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for new proposed method and baselines. If only a subset experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release amonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental sefting/details
Question: Does the paper specify all the training and test details (e.g.. data splits, hyperparameters. how they were chosen, type of optimizer, etc.) necessary understand results?
Answer: [Yes]
Justification: The vital experimental settings are presented the core of the paper, Section And the full details can be referenced in the supplementary document.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, supplemental material.
Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate informa- lion about the statistical significance of the experiments?
Answer: [Yes]
Justification: The experiment results are accompanied by statistical significance tests, i.e., Mann- Whitney test.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to library function, bootstrap, etc.)
The assumptions made should given (c.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
is to report I-sigma error bars, one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CL, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures
symmetric етror bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
Question: For each experiment. does the paper provide sufficient information the computer resources (type of compute workers, memory, time execution) needed reproduce the experiments?
Answer: [Yes]
Justification:
We record the information on computer resources the supplementary document.
Guidelines:
16

---

The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Fthies
Answer:
[Yes]
Justification:
The authors make sure to preserve anonymity.
Guidelines:
The answer means that authors have not reviewed the NeurlPS Code of Fthics.
the authors answer No, they should explain special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g., if there is a special considcration duc to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
This paper is foundational research and not tied to particular applications, let alonc deployments.
Justification:
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has societal impact why the paper does not address societal impact.
Examples of negative socictal impacts include potential malicious unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (c.g., deploy- ment technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimatc point out that improvement in the quality of generative models could used to generate deepfakes for disinformation. On the other hand, it not needed point out a generic algorithm optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional unintentional) misuse the technology.
there are negative socictal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release models, providing defenses in addition to attacks, mechanisms monitor- ing misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g.. pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This paper poses no such risks
Guidelines:
The answer NA means that the paper poses such risks.
17

---

Released models that have high risk for misuse or dual-use should be released with necessary safeguards allow for controlled of the model, for example by requiring that adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification: The creators or original owners of assets (c.g., code, data, models), used in the paper, are properly eredited and are the license and terms of use explicitly mentioned and properly respected.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include a URL.
The name of the license (e.g.. CC-BY 4.0) should be included for each asset. For scraped data from particular source (e.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms of use the package should be provided. For popular datasets, ppaperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide help determine the license of a dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out to the asset's ' creators.
13. New assets
Question: Are introduced in the paper well documented and is the documentation provided alongside assets?
Answer: [NA]
Justification: This paper does not release new assets
Guidelines:
The answer means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their sub- missions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects
Guidelines:
answer NA means that paper does not involve crowdsoureing nor research with human subjects.
Including this information in supplemental material is fine, but main contribution of paper involves human subjects, then as much detail as possible should he included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
18

---

15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: The paper does not involve crowdsourcing research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere to the NeurlPS Code of Ethics the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applica- ble), such as the institution conducting the review,
Declaration LLM usage
Does the paper describe the usage of LLMSs if it an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
LLM not involved other than writing, editing, or formatting purposes this paper.
Guidelines:
answer NA means that core method development this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy http:p////aaaiiaaecc///////rerrs for what should or should not be described.
19

---

A Additional related works
A.1 Evolution strategies
ES 20 23 21 is a prominent class of second-order black-box optimization methods, Compared to zeroth- methods that leverage probability distribution models in their simplest forms and usually amount to finite- difference approaches [52]. ES is considered second-order because they are closely related to the quasi-Newton method [53] they typically learn covariance matrix that converges the inverse of the objective's Hessian.
Wierstra al. 20 proposed the Natural Evolution Stralegy (NES), which amounts a stochastic natural gradient ascent. Since its publication, several NES variants have emerged [54] 28 -3 addressing the issue with the Fisher information matrix using different parameterization strategies. Both Akimoto et al. [40] and Glasmachers et al. [28] derived the rank-jI update of CMA -SEESSAAAA from natural gradient, bridging gradient search and state-of-the-art ES. The primary difference hetween CMA-ES and NES, which makes the former superior, is that CMA-FS leverages two evolution paths that act as the momentum term in gradient search 55 one for updating the covariance matrix, while the other for tuning the isotropic mutation strength.
The covariance matrix the Gaussian distribution is one of the most contributing parameters ES. Through extensive experiments, Hansen et al. (24) showed that with the help of an adaptive covariance matrix, CMA-ES is not affected by coordinate system rotations, and is more resistant to the increase of the condition number. Beyer (26) proved that the covariance matrix gets asymptotically proportional to the inverse of the objective's Hessian, which is preferable as it counteracts the degeneration tendencics of the scarch distribution into a subspace.
The @(n²) space and time complexities regarding the covariance matrix hinder the application of ES to large- scale problems. Several methods have been proposed to remedy the issue. These methods agree on simplifying the covariance matrix, while differing in the details of designing and updating the matrix. Sun et al. [ proposed RI-NES, where rank-one model adopted as the covariance matrix, and the update is completely derived from natural gradient ascent, Akimoto et al. (2 proposed VD-CMA, where covariance matrix is parameterized by the multiplication of a rank-on one model, and a diagonal matrix that captures the scaling of each dimension. They derived the update from natural gradient and, in the meantime, incorporated the cumulation techniques of evolution paths from CMA-ES. The experiments demonstrated that VD-CM& significantly outperformed R RI-NES. The multiplicative model in VD-CMA further pursued in dd-CMA (6), recent method that can well adapt to both separable and non-separable scenarios. straightforward extension to the rank-one model is to apply more than one component to capture multiple promising search directions simultancously, which yields
much more efficient algorithms [7 30 31
Although some recently proposed ESs are versatile across various tasks 23 it may be more appropriate to classify them zeroth-order methods, considering that they mainly focus reducing the variance of the finite difference without utilizing second-order information. On the other hand, while some other methods not explicitly identify as ESs [58 59), they are to some extent inspired by ESs and share many common principles. Notably, Lyu and Tsang (58) proposed a black-box optimization algorithm called Implicit Natural Gradient Optimization (INGO). To avoid explicitly inverting the Fisher, they formulated a surrogate objective where the KL-divergence is attached a penalty term and can he explicitly derived. The authors showed that the gradient with respect to the expectation of the sufficient statistics amounts to the natural gradient of the natural parameter of an exponential-famee distribution. Since the Hessian of an objective primarily characterizes the dimensional dependencies, we focus exclusively on second-orde methods.
A.2 Black-box multi-objective optimization
Multi-ohjective optimization (MOO) has consistently garnered significant attention 34 In black-hox MOO, multi-objective evolutionary algorithms (MOEAs) remain among the most dominant methodologies [ 8 60). This is because MOEAs excel in maintaining a large population of solutions, and can thus approach the underlying Pareto Front (PF) from all directions. Recently, dedicated efforts have been made to extend MOEAs to high-dimensional, or in other words, large-scale MOPs. For instance, IM- IM-MOEA/D (36) integrates the canonical MOEA/D 5 with an inverse model (48) which directly samples in the objective space using Gaussian process models, and then maps into the decision space. Focusing on the complex landscapes in high-dimensional MOPs, LMO-CSO [55] proposes addressing the dimensional dependencies using an efficient particle swarm optimization (PSO) algorithm. Among the three mainstream MOEAs, those based on Pareto dominance and performance indicators fall short of aligning with decision-makers' preferences and scale poorly with the number of solutions [60]. due to the high computational overhead associated with non non-dominated sorting and hypervolume contribution. Meanwhile, decomposition based MOEAs 5 struggle with only a limited number of pre- eeeeealabisee preferences. generalize poorly to the entire Pareto manifold, and encounter difficulties in sharing information across subproblems.
20

---

B More details about Neural-ES
This section is arranged as follows:
Section B.
shows that random samples generated by are in-distribution of T(:: where sp
{
Section B.
derives the vanilla gradient V J.. (0) in (9. and analyzes its major drawback, that is, leads to divergent behaviors of the search distribution.
potentia
Section B.3
justifies the use of the Fisher from the perspective of natural gradient, and derives 10
Section
B.4
further transforms the gradient 100 to the practical update rule (1T) used in the algorithm.
Sections and providd more details on the weight values and the set model.
Section B.
derives the space and time complexity Neural-ES, demonstrating time efficiency.
B.1 Proof of Equation
Proposition A solutoon generated by (8) follows the distribution 7 (::4)), where yp = spo(ex).
Proof. Consider the solution (æ1 ... Zp) where me 0e( 1)(2 ree for all [p]. z (O,,)) I). r (O, 1). To prove x (sp) equivalent proving that (æ1, Lp) follows the joint distribution I No (au, T1, x1-1) where zu a.
Note that for an arbitrary E [p), the expectation of æ; given by:
B[a:;] mo + 00 (x0:s--
[[=] + E(r) vo (Torr-1) mm (atos-i)))
and the covariance matrix is given by:
(a: ma ma(x mo(
= =
R[zz") no [ (1o0 + B6 no (0
or (won--1))) vo (00
Co(wau-1).
which is because E|rz] E[r] 0. as and z are independent from each other. Therefore, for all E [o].
æ;; ~ (<n (1): (o and N(mo (woss Co (1) = No #1, 2;:1)
By the definition of joint distribution, we have:
(#1, 22) (o a1))
(#) 22, ass) No(: (o æ1 Noo a1 x2);
(xi, 22 Fp)
N W- Ii-)).
By this means, show that the reparameterization trick w. ma (200-11 + ar (wros1))( reg (200:-1)) correctly generate in-distribution samples æ T(p)).
B.2
Issue with 9 divergent gaussian mean
Deriving Recall that presents the vanilla gradient of [a( 1x, (f(x)) w.r.t 0, which expressed as:
1 2lnN((2 (p1)
V V.. (0)
Σ E Em((ssp))
f(())
30
dep:
1=1
21

---

where yp; = (a( We are going to derive the term
Wlopee Dep
Notice that the partial derivative of
((;;; sp) = -Il (2: sp.) w.r.f. &p can be derived as:
θ Iu (2; sp)
1
Ón(a; vp)
a (ln dei
(x(w;:) (m: ())
P
0ep т (z: up)
dep
de
Dep:
(12)
Using the above, can derive the partial derivatives 03 /0/p
0J - f(x))
() dx f(z))((2;))) sp)
Σ
Ö In (ansspp) dx
Σ En(mype)
f(x)
a In WinN(at Nax sp.)
ēep:
dep
aei:
ae::
(13)
Eventually, the gradient w.r.t 0 can be derived by directly applying the chain rule, which recovers (9)
in 13) as the vanilla gradient of the log- likelihood of M(x;;; (..) The biggest issue is perhaps that it potentially hinders the convergence of m. which explains why vanilla gradient is not popular and 2) some so-called ES variants proposed recently [23 [2] either update the distribution mean only or update the covariance matrix using gradient- free strategies. For an arbitrary N M(z;; sp) where oo¹ = (m, vec(())") and vec - is the vectorization operator that stacks all the columns of an input matrix [61]. With slight abuse of symbols, we consider in this and the next subsections that the mean m. and covariance matrix C, are constant parameters. The vanilla gradient regarding a random sample æ; can derived as [40]:
The issue. Regarding a random sample æ; E H. we refer to θ In. Nax;
pp:)0pp. (a)::))()))))))))
a In In. sp.) a Cc.)
in N(mssps)) Brn; ((
= (# vec (c.)(a.--))), ";:
Cz: mi) (.)
(14)
Oep:
By eigenvalue decomposition, C, QAAQ", where A is diagonal matrix whose main diagonal consists of the square roots of C, s eigenvalucs, and Q is an orthonormal matrix satisfying Q" =Q æ, can thus be reparameterized by æt; = m. + QAz; where x; ~ NN I). Consequently, we can rewrite
Q
A Q
m, + QAzi -m.) Q A AAQ
QA - QA''Z. 2..
8m:
(15)
whose variance is:
( (Ons((as:)))) =E[(z. Om:
Var
B[z.]
Q a
(z. E[z.)) =
E TT(CT),
(16)
where 2; (Z1y...,77)) 7a)" 14 is the i-th largest eigenvalue of C., =|=']= for all j E [d] is due to the mean of the chi-squared distribution, and Tr(-) denotes the trace a matrix. Intuitively, during the search process, the entropy of the search distribution, ie.. the (log-)determinne Jdeterminant of C. usually decreases. which is likely to increase
Tr(C J. and therefore prevents m. from staying in a stationary point. The phenomenon has been observed
in 20 and rigorously discussed in 26
B.3 Proof of Lemma I& & interpretation of 10
Define x (2; sp) = [[._ (2x spi) and yp (m!" vec C.TYT. Note that treat m: and C, constant parameters instead of mappings from 00:-11 We introduce the definition of natural gradient and its geometric interpretation:
Definition 4 (Natural gradient (32). Consider a distribution (p: that is differentiable worl its parameters sp, with the p.d.f. T(:; (p). The natural gradient of any diffcrentiable function h((p) w.F.t. yp given by
Vh((p)
F(1) Vh((p),
(17)
where and F(()) = (x(gp)))) V In (x; V In क (2; (p | is referred to as the Fisher information matrix w.t. 4p.
Proposition 2 (Geometric interpretation of natural gradient (5)). The natural gradient Vh((p) points in the direction 5 that yields the greatest improvement of h. for given distance between π ((3)p) and T(p 5) in
Kellback-Leibles divergence. More precisely, let up E d be point where [[Vh(())]| 0, then
8" Vh((p) lim |[b((p]]]
2
argmax hlpp + 8). s. (N(((()))) N[ 8))) < 2 (18)
22

---

Of the above, definition 4 provides the form of a natural gradient, and CropositionnSsugeeeeet that natural gradient identities the steepest ascent direction on Riemannian manifolds, and in the meanwhile, guarantees steady changes of m.
We prove Lemma as follows:
Proof. First, the gradient of In (ax; sp) is:
① Iu (w:
④ [(an111))))
Vln.(2 T
V ( (p)
de1
dep
According to the definition of the Fisher,
Wianeee))) ((m1:1))
in Miasspn (p:)) 0 (on,
dep
011
Osp.
F(())= Bey (smup)
(19)
Øln (amnee) a (impee))
I (mppip))) 0 Iln (appopp)))
84p
Dipp
@ln N(m, ((2,
Denoted by Fiz = En(x:)))
the block entry of F(sp) on the i-th row and j-th
dips
column in where $ j. then
Iln Na: (p:) 0 N 4p))
Fg <a;
dw dep)
Oei
P
(( N (I) In Nu da
dep:
Ö4p)
θ In (()(())))))))
Il (ans(i))))))) de
II
N(x;;; ;p;)
04.
04p
(N(e) Deg
(20)
ON @N(xi spd) (N(r)is))))
((zx)
da
dxp
Oep:
04p
01
Bp L.
/ M(
à dip
Na up: dp:
Many sps )an,
Oep,
mj
=Vi=0
=01==0"
which means all the off oif-diagonal block entries are zero, and F(()) turns out to be block-diagonal matrix. 1 Nimuso)) 8 ((eavee) T for all [p), which implics:
Furthermore, notice that F((p.) Em(myy))
Dep
F((P1)
0,
vin = "mo
V In (
0.
0.
0.
F(()))
-1
(21)
F((p,) 8 In Nax:spp
dep
As æ follows T(::p). it unbiased.
Interpretion of 10
We can plug
21 into
10
and apply Monte-Carlo approximation to (J) which
yields:
J.. (0) (x((p)
[(x))
v (In)((:)))
- 2 (fIn)(a):v")))
(22)
80
where sp' (up "p) TT and yis sp ( ). Intuitively, the right-hand side of the above works as follows: first, sample bunch of distributions Tr(p") ), .... Typ?); then, for each of these distributions, sample
23

---

random solution x² ~ T(p)): ): in the meanwhile, O satisfies that pij- p eventually, estimate the natural gradient of In x (xr; w.t.. yp using x². In brief, {10) amounts to aggregating the natural gradient of set of sampled distributions o ((p')..... T(pp), and then proppgaagggg the gradient to 0 through the chain rule.
B.4 Deriving the update rule
22 cannot be directly used in implementing the algorithm, it 1) involves computing F(sp.)". which is intfficient, and does not specify how the gradient back- nock-propagated to 0 through the chain rule. To this end, we first introduce Equation (24) of Akimoto al. [40), which avoids the computation of F(P,))
Proposition 3 (Natural gradient of a Gaussian). Consider a Gaussian random sample æ ~ N((p) where
sp (m,", (e) y², the natural gradient of the log liikeiihood (xx sp.)
F(( Iu [NN((
vee ((x: - -m:)(a - -m.)² -c:)
m. 2x
(23)
depi
Define m., mg (o 1). C. = (I + dada ), or = oo a0ooa-11). vs = vo (2 1), and y: x; m.
According to the chain rule:
=(0m./00)
0 00 F((p:) -1 0 (In.(an;
(or mm)) + Ivec " vee ((oT mi) m.)² C.)
J {m.))²
Oep:
00
=mmm)) yi
(0;
0 vec C₁ J [1:))
O Ce
(v.D? - -c.)
00
do,
Oos
=(Evec Ci/ 00;7
(24)
where J : 12d - Rdx is the Jacobian matrix. In the following, we will expand the above terms one by one:
vec C, vec (n.. C:) vec(I + v.u:)" vec(y.e C:)
20: 0o:
=Th (IT Da D - Tr ((I + v.vi) (T
(25)
TT (y.vr" + Tr y of Tr (1I+
|1 + (oTT.)) od ||0:11²);
and
i vec C:
vec y.y! cc.)
Sv:
=vo(
of 8 m 4
yi 80 y
vec C:)
(26)
of u )
o. (100 1
of (I U. + v
(vec I) + 1z @ v:)
= (o ,4
|1
20 (1 ||: 12)
where denotes the Kronecker product
61 Plugging
and 26 into
24
yields:
(mm))" T yi
+ 20, (Iwvill + (v. y1)² - of 2|| 2 + [|u:11())
do:
(27)
00
+ ((.))
(v. y.u: + |1 - 20 + ||0:11 v)
Consider a sample x², define m.j mo 2 oe
co( ( 1)) fij no ( (mr
and vı = ax. maj-
Eventually, for each j [A] substituting 27
22
yields
B.5 More on the weight values
The rank dependent weights are also known as the use of utility [26225 25 or fitness shaping 220. The rationale
behind the weights in this paper, ie.. ... wx and wi ln( { Inj. is that they decrease
24

---

sub-linearly with the ranks: solutions before the median will be assigned positive weights. while those after (including the median itself) will be assigned negative weights, and no will be assigned zero weight. Separating the weights into positive and negative values is equivalent to using a reward baseline in many practical reinforcement learning methods. Note that in most CMA-ES variants (27, 29 [56,, the second half of { V is discarded by default, especially when updating the Gaussian mean. This is because using inferior solutions implicitly imposes too strong regularity assumption on a black-box objective, which may not be true. A counterexample xNES, which uses the inferior solutions but normalizes their weights to 11/A. However, we adopt all solutions to estimate due to the well-known data-hungry nature of neural networks.
B.6 Design of the set model
We use a feed-forward neural network the model. There are four hidden layers, each with size of 1024. order imput #o into the model, we construct
(n- (4 nila (t-1) 1s
1p-i+1)d as
nil, nil
O,.... 0 20)
(28)
where nil is a placeholder, 1 on the (n (n+)) j) th (1 S j S 12) dimension indicates that the j-th input is not nil, and vice versa. Consequently, the size of the input layer is 2n++ M. The output layer, which outputs (m.." or " size 2d + 1. We use the leaky ReLU function as the activation in the hidden layers, which is defined as ("(w) max(w, 0) + 0.01 min(0, w). In order to generate m,, v; E IR no activation is used in the output layer. To guarantee dr > O, we take the absolute value of the (d + 1)-th output dimension. The set model is trained using Adam, with learning rate of r 10-3. Compared to SGD (stochastic gradient ascent), the adaptive momentum of Adam accumulates historical gradient information, enabling more accurate gradient estimation and significantly reducing the number of samples required.
B.7 Overall computational complexity
Consider set model, i.e., a feed-forwart neural network, with hidden layers, and the sizes of the input, output, and hidden layers are 2n + M, 2d + 1, and 8, respectively, One can derive the space complexity of the network, as well as the time complexities of the feed-forward and brok-proong processes, which are all (lss² + (n + M + d)s). Since we are only interested in and d, we can further simplify the complexity as (n d). We can easily imply that the space complexity of Neural-ESS is O(n + d) well, as the entire search distribution compressed into the model.
The time complexity of Neural-ES stems from two sources: 1) sampling and 2) updating O. To sample a component æ;, it takes (n + d) time to generate a sp. =Poozuu.-)). which corresponds to the feed forward process, and (d) time generate an æ; according 8 To generate complete
(#f ..... 2x" )² S needs
be performed p times. Therefore, overall time complexity for sampling a a ~ a is (n + d) + p = (pm). as n = pd.
To update 0, we need to compute 27
Within the equation, the computation of J Jm.)) de, and J(v,) takes total of (n + d) time, which corresponds to the proppaaatnttiio process; every matrix-vector multiplication
takes (n dd) (dd) ( dd)) time; in addition, the computation of ||y 12 v." y, and || takes (d) time. To up. the time complexity of (27) is O(n d) + O(nd + d²) ((d) O(nd + d²). In addition, (27) needs to be performed p times regarding a complete = (œ T ..... x 1TT which takes
(ne + )²) O(n² +nd)= O(n²) time. Note that ( also represents the time complexity in the update
Bi T as achieves of EPSL 13) @(x²)) is indeed a near-optimun one can expect when updating 0, because
00 Oip
(n²) time only when both yp and 0 scale linearly with n.
A typical counterexample is the application of other existing ES variants, such as xNES, to solve the black-box PSL problem. In this case, the space complexity will increase to O(n²). the time complexity of sampling will increase O(r²), the time complexity of computing the partial derivative (...())) /Oup will increase (n²)) ). and eventually computing the gradient (.) as well as updating O will cost (( time or higher. In contrast, Neural-ES ES achieves ((n) space complexity ((p)) time complexity of sampling a solution. and ((n²) time complexity of updating Therefore, we can conclude that Neural- ES a time-efficient approuch.
C More details on the experiments
C.1 Proposed BBO-PSL benchmark suite
A test problem for assessing the PSL algorithms should satisfy the following two requirements:
25

---

Non-separability, which means argming fa(æ)) o (argmin,, [((xi,..)...................._ arrgiiiin Si(..., ))² for all E [m). This a conditio sine qua non, because separable problems, which are rarely the practical case, cannot assess the ability to handle complex dimension dependency.
Analytical PFs and complex PSs: where the former means that the ground-truth PF and PS should be expressed in closed forms, and the latter suggests that the Pareto optimal solutions should not be, e.g.. linearly distributed.
Surprisingly, few of the existing benchmark suites have satisfied both requirements. For example, the popular DTLZ problems [62 satisfy neither of them: they are separable, with simple PSs located on the boundary of the search space. Despite their complex PS shapes, the UF + and LZ
66 instances are separable. Though
lack analytical PFs & PSs. The LS-MOP [46] focusing on large-scale multiobjective problems does not feature complex PSs, and most of its instances are separable partially separable. An exception would be the recently proposed ZCAT suite +7 which, however, features non simplex- like PFs and falls beyond the scope of this paper.
providing non-separable instances, the bi-objective BBOB problems
4S
For these reasons, we design a novel benchmark suite based on the following theorem, which also the design principle of the popular UF and LZ problems:
Theorem 1 (Li & Zhang 63
Consider a MOP in (7) where
f((x) = (x1) + 8 (x111 ((((()))
(29)
i I,.., M. x2 = (Xi,.., XM- I11 = (x)),..,Xa). a E X, X 1o,= C R", S
- [os, + 12, lr 11 (oj,
10--M+1, and pote IT [I-MM+1as (aa, by] + RR. Assume 1)
fod reaches the minimal of O, which means 1 (x) 2 pd (0) = O. Va RR, and 2) the following MOP has a PF and PS, which are denoted by PF and PS",
minimize T [o..bl)
(f (.., S (a:))²²
(30)
Then, 1) PF of Fx) PF", and 2) the PS of F(x) is (21,, " f"')),, PS*).
Of the above. 1." f and fdis are referred to as the PF, PS, and distance functions. as they determine the shapes of the PF and PS manifolds well as the difficulty of optimizing F(x)), respectively. With Theorem analytical PFs and complex PSs can be implemented by designing proper s." and f and non-separability can be implemented by introducing non-separable functions that achieve the minimum of zero
O.
We design instances in the proposed benchmark suite as follows. All the instances employ the rotated ellipsoidal function as the distance function, which, considering an input æ11, expressed as:
(T R) = RAR" arr,
(31)
where R is (n -M d--menssional orthonormal matrix, and A is (n M + 1)-dimensional diagonal
matrix whose i-th diagonal entry is c Particularly, c called the condition number (3,, which denotes the ratio between the largest and smallest cigenvalues of a Hessian, and typically controls the optimization difficulty; a larger yields more difficult optimization problem. R can be constructed by applying the Schmidt process to random square matrix with standard normally distributed entries. All the bi-objective instances take the following form:
[f(((m) f (((1 - (m)); R) (x1) + lrd (x11- - f (21); R')
F(a)
(32)
where the PF function is fr (211 (21, vz1)" All the tri tbjeotiee instances take the following form:
At (2), + fod (œ11 (r1, R)
H (x1. = cos (0.5wzz) cos(0.5mmz)))
F(x)=
ftr (21, + fol [11 ( x2); R')
where
1 (x1, x22 = co.(0.5 sin(0.5mmr2))
19 (21, + f (œ11- (m, 22); R")
Ard (x), = sin50.)))
(33)
and R. R' and R" denote three different rotation matrices.
Eventually, eight instances, Fl to F8, are listed in Table 2 Fl to F4 are bi-objective problems while F5 F8 are
trieotvecivve. The search space for the bi-objeetive instances is [0, 1] x [-1,1] and is [0, 1J2 x 2- for the tri-objective ones. The condition number c 100 in Fl to F4, and is 2 in F5 to F8. Their PS manifolds in the (2i 22 r²) subspace are displayed in Fig. g In short, the proposed suite features non-separable and anisotropic distance functions, covering diverse nonlinear PS manifolds.
26

---

Table 2: The proposed BBO-PSL benchmark suite.
ID
M
f
(x1 0.75(1 - 1/)a1 cos(22ta ), if i is even.
F1
sin (2n(2)11) otherwise
2
(v/s (= 2x1 + 0.5), if is even,
F2
(x1)
1 -0.5, otherwise
511
F3
(1 = (1 - (2x1 1)3 + (2 - 1)
F4
(x1) 2x 0.3+1.36.
F5
(2) 222 = 0.56; cos( n 2 COS(TT])) I2
F6
ft (21 22) = &1 (2î - 12) + + (1 x9) -
F7
(x1,32) (12((i(((T))))
2) + (2--22)(66) (sin TT 2))
3
F8
22) = -
82 -1)
f. (-) denotes the i-th output of f., v; == and &: = n-i-1
m-
(a) Fl
(b) F2
(c) F3
(d) F4
(c) F5
(f) F6
(g) F7
(h) F8
Figure 8: Ground-truth PS manifolds of the proposed problem instances.
C.2 More about the configurations
Compared algorithms. The two MOEAs used in the experiments, namely LMO-CSO and IM-MOEA/D, are implemented in the open-soure PlatMMO library [64] Their hyperparameterss are configured to their default values. The set models of EPSL and EPSL-RI use the same hidden-layer structure Neural-ES. Berouuli smoothing, antithetic sampling, and smooth Tehebycheff EPSL are deactivated
fair comparison. The most 4 + [3 In n a typical setup most ES algorithms [19 28), reduce half due the large number of sampled subproblems (g(:; a 2 We the popular Penalty-base based Boundary Intersection (PBI) as
dominating hyperparameters Neural-ES are p. N, and X. While
We will discuss p N in Section C.4
the aggregation function:
Dr 1 (u (( all |al
Pht phi
(F; Dr + pD2, where
D2 = F + D.all:
F is a point in the objective space, B reference point satisfying u S (inf (1x) inf fM (x)) Dj is the distance from u to the projection onto a, and Da is distance from a. More specifically, by minimizing D,. one approaches the underlying PF, while minimizing D2, which serves as a penalty, aligns
27

---

with the preference balpha. The parameter p > 0 trades off Dj and Dz. Although the Tchebycheff function [8] is also widely used, we have not considered because we observe that its decrease solution diversity when solving tri-objective MOPs, which further leads to premature convergence of the set model.
Performance metries. The Inverted Generational Distance (IGD) 49 is used for assessing the algorithms'
performance, with a smaller value indicating better performance. Denoted by PF and PF subset of the ground truth PF, and PF an approximation by an algorithm, the IGD of PF is defined as the total distance from every point in PF to PF, and can expressed as:
IGD(PE, PF)
Σ PEEPF
min F
PF
Intuitively, IGD measures the distance from each point F in the ideal set PF. the closest F in the approximation set PF. When imposing ideal properties on PF, the IGD is minimized only when PF spreads widely and evenly. Otherwise, there will always he some F E PF that remains distant from PF, resulting in large IGD.
Additionally, hypervolume (HV) is also used to evaluate algorithms in Section 5.3 where IGD values are unavailable. It measures the volume of a specific region in the objective space, which is between a solution set and a reference point, with a larger value indicating better performance. is given by:
(S) = Vol Ra E S, F(x) - - }
where S solution set, is reference point in the F F-space, Vol computes the volume of a set, and < denotes Pareto dominance. The volume always increases when the solution set approaches the ground truth PF closely or spreads more widely across the PF, as hoth cases stretch the region.
Experiment configurations. We limit the size of PF be 300 and 900 for the bi- and tri-objective instances, respectively. Each algorithm repeatedly solves each BBO-PSL instance 21 times, using a different random seed cach time. For the three PSL algorithms, we evaluate their testing performance. Specifically, we input large number of preferences into a trained set model, which yields the same number of distributions; from each distribution, we randomly sample a solution, the collection of which form the testing solution set; finally, we measure the IGD value of the set and refer it as the testing IGD value. The number input preferences, i... the size of the testing solution set, is 900 and 4950 for the bi- and tri tri-objective instances, respectively. The MOEAs directly optimize the same number of solutions, and the IGD value of the final population is measured. Meanwhile, we track the convergence of the PSL algorithms during training. in terms of the IGD values of the training solution set. The training solution at certain iteration comprises the hest solutions found so far for a set of evenly distributed subproblems. In light the computational efficiency, its capacity is slightly smaller than the testing solution set, which is 100 and 300 for the hi- and tri-ohjective instances, respectively. To examine the scalability of the algorithms with increasing dimensions, we vary n from 32 to 1024. The evaluation budget is 4000n and 2000n. for bi- and tri-objective instances. The bi-objective instances given a larger budget because their condition numbers are higher, which implies greater difficulties.
C.3 More experiment results
Numerical results. The numerical results for Fig. 4 are listed in Table 3 The Mann- WhinneW test at a 95% confidence level is applied to statistically compare the algorithms. Neural-ES performs the best in all the instances, ranging from n = 64 to 1024. The results conform to those in Fig. 4 Note that the results of LMO-CMO and IM-MOEA/D are obtained through directly optimizing the problems, while the results of the three PSL algorithms represent their test-time gencralization outcomes. This indicates that even the testing performance of Neural-ES surpasses that of the MOEAs terms of optimization results.
Approximate PF manifolds. Fig. D demonstrates the approximate PF manifolds by all five algorithms, on all the BBO-PSL instances. In all cases. the PFs approximated by Neural-ES dominate most areas of those approximated by the other algorithms, which conforms to Fig. 5 It appears that F3 and F4 are more challenging than FI and F2, and F8 more demanding than F5 to F7, indicating that the shape of the PS manifolds significantly influences the optimization difficulty. It can be observed that Neural-ES closely approximates the ground-truth PFs of all bi-objective instances, which highlights the learning capability of black-box PSL algorithms, particularly Neural-ES. However, the approximation results are inferior in the bi-objective instances, due their significantly higher condition number. This confirms that the condition number of the distance function notably impacts optimization difficulty, which existing benchmark suites have neglected. Once again, the approximate PFs of Neural- ES and the other two PSL algorithms are produced during testing time, while those of LMO-CSC and IM- MOEA/D are direct optimization results. This suggests that the generalization ability of Neural- ES surpasses that of some state-of the MOEAs.
28

---

DAMPAAT
SOMPAA
DPS
EPN
ENS
EPA
SIUREE PP
(a) Fl
(b) F2
(c) F3
(d) F4
L4000
4
4
2010
4 14000
(e) F5
(0) F6
(g) F7
(h) F8
Figure 9: Testing-time PFs approximated by the algorithms, with n 256. For better visualization, only non-dominated solutions are presented for F5 F8.
PS RPRL
PS EPSI
PS SPSL. SPSL-RI
PS SPSL,
EPSI-RI
SPSI-RI R1
Noural ES
Neural 25
Neural ES
Neural ES
(a) FI
(b) F2
(c) F3
(d) F4
Figure 10: Testing- Testing-time PSs approximated by the PSL algorithms, on F1 to F4, with n = 512.
Approximate PS manifolds. We focus on PSL algorithms' ability to learn manifolds. LMO-CSO and IM-MOEA/D are excluded as they do not possess such ability. The approximate PSs on the bi-objective Fl to F4 are plotted in Fig. [00 Neural-ES generates the most accurate PSs in all four cases, as the red curves align the closest the ground- truth. It can be observed that F3 and F4 are more challenging than Fl and F2, they have higher levels of non-lincarity. The approximate PSs on the tri-objective F5 F8 plotted in Figs. Figs.ITtof 14 Neural-ES and EPSLL-R generate closer and more uniform approximations than EPSL. in four cases. Although Neural-ES only slightly outperforms EPSL,RR in the solution space, their performance gap is significant, because most solutions of Neural-ES dominate those of EPSL-RI in the objective space, as shown in Fig. 9
Convergence. Fig. 15 demonstrates the convergence trajectories of the three PSL algorithms. LMO-CSO
and IM-MOEA/D are not involved due to their differing sample sizes during each itcration. The results conform to Fig. 3 Noticeably, EPSL-RI converges slightly faster than EPSL in Fl to F4; however, instances F5 to F8 do not reflect its advantage. This is because Fl to F4 have much higher condition number of 100, than F5 F8, which have condition number of 2. Evidently, dealing with the dimensional dependencies becomes more challenging as the condition number increases. Once again, existing benchmark suites overlook this issue, necessitating the proposal of the BBO-PSL suite. In all the instances, Neural-ES establishes a performance lead over EPSL and EPSL-RI during the early stages, extending the performance gap throughout the search process, underscoring its superiority in handling dimensional dependencies.
29

---

PS EPSL-RI
PS Neural-ES
EPSL
0.0
(a) F5, EPSL
(b) F5, EPSL-RI
(c) F5, Neural-ES
Figure 11: Testing-timg time PSs approximated by the PSL algorithms, on F5, with n 128.
PS Neural-ES
EPSL
EPSLRI
15
42
4.0
(a) F6, EPSL
(b) F6, EPSL-R1
(c) F6, Neural-ES
Figure 12: Testing-time PSs approximated by the PSL algorithms, on F6, with n 128.
C.4 Sensitivity analyses
Intuitively, the number of partitions, p. and number of sampled subproblems, N, are closely linked to the performance Neural-ES. To determine their optimal setups as well as to observe how they vary with the increasing problem dimension we conducted the following two independent experiments:
Investigating p. apply p (64, 16, 8,4) solve F4 and F6, and meanwhile, fix N [5n/16].
Investigating N: apply N E {10, 40, 80, 160, 320) to solve F4 and F6, and meanwhile, fix p 16.
In the experiments, the dimensions of F4 and F6 range from 64 to 512. The two selected instances serve as the representative bi- and tri-objective 3BO-PSL instances, which are expected to provide the most neutral observations. The experiment results, in terms of the final IGD values, are recorded in Fig. 16
We can make the following observations:
Regarding Figs. [6(a) and [6(b] number of partitions p positively correlates to the problem dimension n. This is reasonable as a larger p yields more decomposed Gaussians, and therefore a stronger ability to handle dimensional dependencies. The results validate the efficacy of the novel dependency handling technique presented in this paper (i.c., Equation ()).
Regarding Figs. [6(c) and [6(c] the highlighted main diagonals of the figures indicate the trade-offss between N and Neural ES: either too large too small sample size worsens the performance. This
is rcasonable because, the one hand, too large N consumes significant portion of the evaluation budget; on the other hand, too small N results inaccurate gradient estimates.
To balance computational efficiency and performance, we adopt = 16. It can be observed from Figs. [6(c) and [6(d) that the optimal N is roughly [5n/16]. To date, we have justified the hyperparameter settings presented in this paper.
30

---

PS Neural-ES
EPSL
ECSL-R
(a) F7, EPSL
(b) F7, EPSL-RI
(c) F7, Neural-ES
Figure 13: Testing-time time PSs approximated by the PSL algorithms, on F7, with n 128.
EPSL
EPSLRI
Neural-ES
-1
12
12
(a) F8, EPSL
(b) F8, EPSL-RI
(c) F8, Neural-E3
Figure 14: Testing-time PSs approximated by the PSL algorithms, on F8, with n = 128.
C.5 Runtime analyses
We have conducted extensive evaluation of computational efficiency across all compared algorithms using 1024-dimensional instances (F1 and F5) from BBO-PSI. suite. The evaluation budget remained the same as hefore, at 4, 000n.
and 2, 000n for F5. The results are reported Table
Comparing PSI. algorithms with traditional MOEAS, we can ohserve that the former run faster, the latter. This is because most MOEAs require advanced techniques to maintain population diversity and address slow con- vergence, such as the RVEA-based en- vironmental selection in LMO-CSO and the Gaussian process modeling in IM- MOEA/D. Most, not all, of these techniques come with high overhead scale poorly with increasing population sizes. In contrast, the PSL, paradigm does not rely on any of these techniques. Despite introducing neural networks. the PSL algorithms turn out run surprisingly faster, or at least equally fast, when compared traditional MOEAs.
least as fast, as
Table 4: Wall-clock time (in minutes).
LMO-
IM-
EPSI.-
Neural-
EPSL
CSO
MOEA/D
ES
F1
16.9
17.4
7.8
12.2
17.4
F5
23.9
15.2
5.1
7.4
9.5
While Neural-EE ES shows marginally higher runtime compared to the other two PSL algorithms. due its more complex search distribution, it achieves much faster convergence terms of number of generations required. The faster convergence likely offsets its marginally higher computational cost per generation, making it advantagcous for practical applications.
31

---

EPSLIT
[P5.311
Neatli
(a)
(h)
(d)
[PA-01
ESLL01
(c)
(
(g)
(h) F8
Figure 15: Convergence curves regarding the logarithm of the IGD values on the training solution sets, with n 256.
3.97
3.54
5.77
NA
6.67
4.59
5.08
6.29
3.92
3.97
4.59
11.
25.2
359
849
6.56
9.2
5.75
11.6
510
128
88
8.42
6.60
5.40
8.85
128
5.23
4.46 4.44
5.09
10.6
238
8.69
4.72
4.60
4.37
5.62
731
256
9.22
4.37
5.00
8.32
11.3
2.25
4.32
4.75
5.72
7.90
512 25g
7.59
5.97
5.45
5.00
5.75
145
256
30
8.27
4.75
6.96
5.7
512
65
3:38
1:46
10.3
14.5
512
3:30
4.43
5.05
8.45
8.71
8.73
3.04
8.
8.07
7.46
7.52
6.13
5.96
6.64
5.27
5.05 4.3
(a) P versus 71, on F4
(b) versus n, on F6
(c) N versus on F4
(d) N versus n, on F6
Figure 16: The acquired IGD values (x10- in sensitivity analyses.
32

---

Table Average training testing IGD values (x10 and standard deviations (x 1x10 in paren- theses).
7
LSO-CMO
IM-MOEA/D
EPSL
EPSL-RI
Neural-ES
2
12.45 (47. 98.9
5.141
7.202 (4.887)
6.946
6.099 (8.028)
26
15.77 (41.37)
5.570 (19.31)
7.187 (4.208)
8.216 (8.306)
4.797 (8.051)
2
20.38 (30.18)
7.168 (21.02)
8. 8.252 (4,220)
9.453 (8.892)
3.997 (6.794)
F1
2
22.11 (22.86)
7.231 (17.35)
9.758 (4.906)
12.10 (7.956)
3.967 (6.890)
2
24.05 (9.820)
8, 8.082 (12:90)
11.04 (3.094)
14.74 (6.150)
3.997 (4.171)
210
26.62 (5.969)
10.69 (9.372)
12.71 (3.838)
16.88 (6.946)
5.336 (9.748)
2
4.993 (17. 411
6.526 (10.84)
7.016 (8.088)
4.588 (12.11)
5.457 (16.48)
26
6.739 (19.71)
7.863 (14.81)
7.449 (7.532)
6.588 (15.80)
3.883 (11.72)
F2
2
8.445 (14.83)
8.726 (12.97)
8.773 (4.078)
7.581 (12.42)
3.004
(3.096)
2
10.54 (12.55)
9.788 (10.07)
11.79 (7.447)
11.95 (17.86)
3.489
10.869)
29
15.13
(21.87)
13.20 (34.25)
13.59 (4.785)
15. (13.01)
3.727
(5.935)
210
17.65 (16.18)
13.20 (15.28)
17.05 (6.887)
19.12 (12.04)
4.889
(5.704)
2
7.497 (12.54)
13.52 (9.551)
9.410 (7.694)
6.629 (7.717)
6.692
(11.20)
26
9.351 (5.085)
14.30
9.486 (7.293)
7.477 (9.134)
5.771 (8.695)
2
11.83 (3.838)
16.12 (4.709)
10.63 (6.437)
8.397 (7.823)
5.986 (9.637)
F3
2
14.24 (6.441)
16.61 (4.133)
11.68 (5.749)
10.13 (8.249)
5.863 (6.778)
29
16.61 (8.691)
16.66 (2.910)
12.32 (5.782)
12.81 (7.730)
6.078 (4.269)
2
19.13 (5.408)
16.76 (1.821)
14.22 (6.943)
15.81 (9.210)
6.654 (6.296)
2
4.647 (39.05)
20.90 (13.96)
15.85 (18.37)
12.23 (13.27)
6.561 (12.06)
2
8,420 (28. (28.84)
21.74 (11.57)
14.44 (14.62)
12.67 (18.09)
5.685 (9.287)
F4
2
17.70
(3)
23.77
21.53 (9.107)
15.71 (23.61)
6.558
28
23.87 (16.02)
24.73 (6.528)
23.66 (7.163)
19.13 (12.07)
7.275 (8.371)
29
27.69 (7.691)
25.25 (4.206)
24.64
(4.299)
22.88 (15.45)
8.728 (8.730)
2
29.15
25.42 (7.191)
26.29
26.65
10.65
(10.97)
2
4.190 (2.761)
12.31 (12.62)
7.233 (6.933)
5.494 (6.135)
3.261
(6.037)
2
5.815 (4.487)
19.31 (15.53)
8.065 (9.451)
6.052 (6.118)
3.090
(5.485)
2
8.869 (7.198)
24.87 (18:05)
7.806 (6.450)
7.234 (7.670)
2.793 (4.118)
F5
28
12.22 (17.98)
32.98 (27.00)
8.830 (6.403)
7.982 (6.154)
3.147
(3.708)
29
14.39 (16.78)
36.49 (21.15)
8.976 (6.300)
9.016 (11.79)
3.370
(4.695)
2
16.79 (21.62)
39.58 (25.25)
9.33 9.383
9.764 (7.775)
3.937
(4.777)
2
4.528
(2.817)
13.23
(13.22)
7.454 (10.10)
5.991 (6.928)
2.752 (3.492)
26
6.351
(5.353)
20.27 (16.66)
8. 8.500 (8.668)
7.368 (6.723)
2.764 (3.685)
F6
27
9.792 (13.22)
33.81
(129.79)
8.930 (5.306)
8.382 (8.728)
2.635
(4.445)
28
14.04 (23.30)
32.74
(18.43)
9.503 (4.951)
9.297 (6.391)
2.917
(3.655)
29
18.06 (24.40)
37.48
(18.
10.10 (5.686)
10.33 (6.853)
3.261
(3.124)
220
20.92 (41.97)
39.94
(11.84)
10.32 (4.074)
11.31 (13.02)
4.041
(4.753)
2
4.633
13.38
7.998 (6.516)
5.918
4.026
(6.118)
26
6.498 (3.846)
18.77 (16.76)
8. 1.0 (7.429)
6.694 (4.773)
4.166
(3.573)
F7
2
9.592 (13.00)
28.54 (50.67)
8.581 (4.287)
7.147 (6.947)
4.213
(3.333)
28
12.69 (17.68)
31.79 (17.96)
9.277 (7.330)
8.386 (7.677)
4.330 (2.266)
29
17.31
(27.64)
34.64 (12.73)
9.585 (6.458)
9.084 (7.185)
4.523 (2.566)
210
21.01
(23.59) 59)
37.25 (17.04)
9.726 192)
9.935 (8.296)
4.833 (1.760)
2
4.833 (3.210)
12.41 (22.91)
7.790 (7.609)
6.726 (7.825)
5.105 (5.517)
26
6.765 (5.986)
18.51 (12.94)
8.413 (9.989)
6.936 (6.334)
4.472 (4.868)
2
9.127 (8.242)
25.61 (22.68)
8.663 (7.097)
7.340 (0.793)
4.220 (3.094)
F8
28
12.44 (14.13)
30.38 (15.34)
9.471 (7.997)
8.243 (6.328)
4.347 (3.694)
29
16.69 (22.96)
34.29 (16.27)
9.076 (5.806)
8.655 (6.327)
4.477 (2.702)
210
19.06 (42.91)
36.74 (17.06)
9.462
(4.698)
9.316 (5.610)
4.939
(3.172)
33


Reviews:
[{"review_id": "kPAk4HeZ7E", "reviewer": "Reviewer_zJXa", "summary": "This paper proposes a NES-based optimization method to solve high-dimensional black-box multi-objective optimization problems (MOPs) via Pareto Set Learning (PSL). PSL maps user preferences to Pareto-optimal solutions using a set model (a neural network). Neural-ES partitions high-dimensional variables into segments, models dependencies via conditional Gaussian distributions parameterized by a neural network, and uses a novel natural gradient estimator for training. Notably, the method scales PSL to over 1000 dimensions.", "questions": "Please refer to 'weaknesses' in Strengths And Weaknesses.", "limitations": "Yes.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. The proposed method integrates ES with neural networks to capture complex dependencies in high-dimensional spaces, addressing a critical gap in black-box PSL.\n\n2. The proposed method demonstrated scalability to 1024-Dim problems, which significantly beyond prior PSL methods.\n\n3. The authors leverage natural gradient estimation and Fisher information for stable convergence, with proofs in supplementary materials.\n\nWeaknesses\n\n1. Neural-ES requires training a neural network and running ES evaluations. The wall-clock time vs. traditional ES/MOEAs is unreported, raising concerns about practicality.\n\n2. Complexity claims of the proposed method focus on asymptotics but ignore constants (e.g., network inference costs).\n\n3. BBO-PSL functions are synthetic. Real-world high-dimensional MOPs (e.g., neural architecture search) remain unvalidated.\n\n4. The method assumes unconstrained MOPs, but real-world problems often include constraints. The limitation is acknowledged but unaddressed.", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate the reviewer's thoughtful and insightful feedback. Below, we provide detailed responses to the key concerns raised, with new experimental validations incorporated in the revised manuscript. \n\n------\n\n> ### **W1:** Runtime assessment\n\n**A1:** We have conducted an extensive evaluation of computational efficiency across all compared algorithms using two 1024-dimensional instances (F1 and F5) from our BBO-PSL suite. The evaluation budget was the same as mentioned in Sec. C2 of the supplementary document, which is $4000n$ for F1, and $2000n$ for F5. The runtime results are as follows: \n| Algorithm | F1 ($n=1024$) | F5 ($n=1024$) |\n|:------------:|:------------------:|:------------------:|\n| LMO-CSO | 1019s (16.98min) | 1432s (23.87min) |\n| IM-MOEA/D | 1044s (17.40min) | 911s (15.18min) |\n| EPSL  | 486s (7.76min) | 308s (5.13min) |\n| EPSL-R1 | 731s (12.18min) | 443s (7.38min) |\n| Neural-ES | 1042s (17.37min) | 567s (9.45min) |\n\nWe make the following observations. First, the PSL paradigm demonstrates remarkable computational efficiency. Both EPSL and EPSL-R1 achieve significantly faster runtimes (2-3×) compared to state-of-the-art MOEAs. This is because most MOEAs require advanced techniques to maintain population diversity and to deal with slow convergence, such as the RVEA-based environmental selection in LMO-CSO and the Gaussian process modeling in IM-MOEA/D. Most, if not all, of these techniques come with high overhead and scale poorly with increasing population sizes. In contrast, the PSL paradigm does not rely on any of these techniques. Despite employing neural networks, the PSL algorithms turn out to run surprisingly faster, or at least equally fast, when compared to traditional MOEAs. \n\nWhile Neural-ES shows marginally higher runtime compared to the other two PSL algorithms, due to its more sophisticated distribution modeling, it achieves dramatically faster convergence in terms of generations required. As shown in Fig. 15 of the supplementary document, Neural-ES needs only one-fifth to one-tenth the number of generations to reach the same IGD performance level as EPSL and EPSL-R1. This faster convergence likely offsets its marginally higher computational cost per generation, making it advantageous for practical applications.\n\n------\n\n> ### **W2:** Complexity analyses\n\n**A2:** We analyzed the overall computational complexity, including the network inference costs, in Sec. B-7 of the supplementary document. In brief, Neural-ES achieves $\\mathcal{O}(n)$ space complexity and $\\mathcal{O}(n)$ time complexity when ignoring the network inference cost, as well as $\\mathcal{O}(n^2)$ time complexity when taking into account the cost of the network. This complexity is identical to EPSL and EPSL-R1. This matches the theoretical lower bound for neural network-based gradient estimation, where the chain rule operation $\\frac{\\partial}{\\partial\\theta}\\left(\\frac{\\partial\\mathcal{J}}{\\partial\\varphi}\\right)$ inherently requires $\\mathcal{O}(n^2)$ computations when the distribution parameters $\\varphi$ and the network parameters $\\theta$ both scale linearly with dimension $n$. In light of the superior performance of Neural-ES, we believe the benefits justify the costs. \n\n------\n\n> ### **W3:** Real-world problems\n\n**A3:** In Section 5.3, we presented an experiment on a real-world trajectory planning problem. However, in response to the valuable suggestion to validate our method on more practical problems, we have supplemented our work with new experiments focusing on multi-objective unmanned aerial vehicle (MO-UAV) navigation, which is featured in the Meta-BBO-v2 library [1,2]. This problem consists of three competing objectives, which are to be minimized: \n\n* the distance from the UAV's final location to the destination, penalized by the trajectory’s non-smoothness; \n* the risk of collision with threats, such as ground-to-air missiles; \n* an altitude penalty to maintain safe and cost-efficient flight. \n\nThese objectives are thoroughly elaborated in [3]. By varying the control parameters of the terrain, we have developed three easy, moderate, and hard instances, which require a UAV trajectory of 10, 20, and 40 waypoints, respectively. Because each point comprises x, y, and z coordinates, the overall problem dimensions are 30, 60, and 120, respectively. MO-UAV serves as a typical real-world scenario where complex dimensional dependencies arise from temporal correlations between waypoints, which makes the navigation highly challenging. \n\nAfter 21 independent runs, we report the brief results in terms of the average HV values (higher is better) and their standard deviations as follows: \n| Instance | Dimensions | Neural-ES (HV) | EPSL (HV) | EPSL-R1 (HV) |\n|:---------:|:----------:|:--------------:|:------------:|:-------------:|\n| Easy | 30 | **3.397±0.279**  | 2.471±0.356 | 1.577±0.149 |\n| Moderate | 60 | **3.443±0.452**  | 2.399±0.372 | 1.520±0.267 |\n| Hard | 120 | **2.442±1.143**  | 0.716±0.197 | 0.842±0.139 |\n\nThe significantly higher HV values across all difficulty levels confirm Neural-ES's effectiveness in practical multi-objective optimization scenarios.\n\n------\n\n> ### **W4:** Limitations\n\n**A4:** We fully acknowledge the importance of constraint handling in multi-objective optimization, as rightly emphasized by the reviewer. Constraints are indeed a fundamental challenge in real-world applications. However, constraints pose particularly unique difficulties for Neural-ES, which stem from two core algorithmic characteristics: \n* __Manifold discontinuity__: Constrained MOPs often produce fragmented Pareto fronts with isolated regions. This disrupts Neural-ES’s key assumption of a full mapping from preference vectors to Pareto optimal solutions.\n* __Sampling breakdown__: The sampling in Neural-ES, where a segment $x_i$ depends on preceding $x_0,...,x_{i-1}$ becomes unrecoverable if early segments frequently violate constraints; besides, gradient-based training also fails when constraints block feasible ascent directions. \n\nWhile we acknowledge this critical gap, addressing it thoroughly would require further algorithmic changes, which may be too complex to cover in detail within the scope of this brief paper. We genuinely appreciate this critique and will prioritize it in subsequent research.\n\n------\n\n> ### **References**\n\n[1] Z. Ma, et al. \"Metabox: A benchmark platform for meta-black-box optimization with reinforcement learning.\" NeurIPS 36 (2023): 10775-10795.\n\n[2] Z. Ma, et al. \"MetaBox-v2: A Unified Benchmark Platform for Meta-Black-Box Optimization.\" arXiv preprint arXiv:2505.17745 (2025).\n\n[3] MA Shehadeh, J Kůdela. \"Benchmarking global optimization techniques for unmanned aerial vehicle path planning.\" Expert Systems with Applications (2025): 128645."}, {"review_id": "7xQLsGXn63", "reviewer": "Reviewer_37Ck", "summary": "The paper proposes a new method for black-box multiobjective optimization based on Pareto Set Learning (PSL). The approach relies on a Evolution Strategy (ES) scheme to conduct the optimization for specific objectives preference (trade-off), which search for specific preference is fed back to PSL for updating the general PSL set model. The PSL set model provides the sampling distribution used by the ES for the given preference. Results show significantly better performances that several other multi-objective optimization techniques applied to black box problem, including approaches based on PSL.", "questions": "- Can you provide a more detailled assessment of the overall efficiency of the proposed Neural-ES approach regarding the number of samples to evaluate for converging to good set compared to other approaches?\n- I think there is a notion of coverage over the Pareto front is important for the current problem, but this is not much explored, metrics used are IGD and hypervolume. Can you explain how these metrics relate to the fact that the Pareto front is well covered (or not) from a given optimization run?", "limitations": "There is not much direct impact on society coming from the current work, so no discussion is required here. The Limitation paragraph in the Conclusion (sec. 6) consider constraint in multi-objective optimization as a future work is quite out-of-scope and seems to be a general problem with PSL. I am not sure it is much relevant to mention this in the paper.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n- A quite straightforward and sound proposal to extend PSL to black box multi-objective optimization.\n- Detailed presentation of the method, with enough details that should allow to reproduce the method from the paper rather easily.\n- The paper is quite well-written and easy to follow.\n\nWeaknesses:\n- Black box optimization usually requires more potential solutions (samples) to be evaluated for achieving good enough solutions compared to cases where analytical equations of the objectives are available. Except for Fig. 3, which shows convergence over generations in a specific case, this is not considered in much details in the paper, while being an important aspect to consider, in my opinion.\n- The experiments are not always easy to follow and the figures presented are not that obvious to understand. For instance, from Fig. 6 (a)-(c), I am not sure about conclusions to have.\n- Likewise, the practical application in Sec. 5.3 is not super convincing, the results in Fig. 7 showing simply that the proposed Neural-ES approach is skewed toward low f2 values, while EPSL-R1 is the opposite. Many of the results presentated appear more visual/qualitative than anything.", "was_revised": true, "final_justification": "I would like to thanks the authors for their answers and extra experiments conducted. If not already done, I am advising strongly to update the paper according to the discussions with the reviewers. I am satisfied with the answers and I am maintaining my evaluation.", "rebuttal": "Thank you very much for your supportive recommendation and insightful comments. We are thrilled that you found our work straightforward, sound, and easy to follow. We will include all the suggested details and explanations in the revised manuscript. Below, please see our detailed responses.\n\n------\n\n> ### **W1: & Q1**: Number of samples required for convergence.\n\n**A1:** The total number of samples is proportional to the total number of evaluations. In each generation of the Neural-ES, EPSL-R1, and EPSL, the number of solutions drawn is $N\\lambda$, where $N$ and $\\lambda$ refer to the number of sampled preferences and the number of solutions for each preference, respectively. Thus, the total samples across all generations are $N\\lambda\\times$ (total generations). As specified in Sec. C.2 of the supplementary document, $N=\\lfloor5n/16\\rfloor$ and $\\lambda=2+⌊1.5\\ln⁡{n}⌋$, and the total number of solutions is $4000n$ and $2000n$ for the bi- and tri-objective BBO-PSL instances, respectively. \n\nWe plotted the convergence curves for all eight BBO-PSL instances in Fig. 15 of the supplementary document. To achieve the same IGD performance, Neural-ES always uses much fewer generations (and thus fewer samples) than EPSL and EPSL-R1. This highlights its sample efficiency, critical for black-box optimization.\n\n**Role of Adam Optimizer** Further efficiency gains arise from the use of Adam for training the set models in all three algorithms. Adam’s momentum-based gradient accumulation improves gradient estimation accuracy over time, which reduces the need for excessive sampling.\n\n------\n\n> ### **W2**: Observations from Fig. 6\n\n**A2:** We sincerely apologize for any confusion caused by the insufficient explanation of Figure 6 due to space constraints. Figs. 6(a) to (c) depict the approximate Pareto sets in the first three dimensions of the search space (i.e., the $x_1-x_2-x_3$ subspace of the $x$-space). In these figures, the gray surface is the ground-truth PS, while the colored points indicate the sampled solutions by different algorithms during the testing time. The solutions produced by EPSL appear scattered and poorly aligned with preference vectors, while both EPSL-R1 and Neural-ES demonstrate better organization. \n\nFig. 6(d) depicts the ground-truth Pareto front (i.e., surface with gray contours) and the algorithm approximations in the objective space (i.e., the $F$-space). For clarity, we excluded dominated solutions from visualization. This provides a more rigorous performance assessment because it directly evaluates solution quality in the objective space and reveals dominance relationships not apparent in decision space. While Neural-ES and EPSL-R1 appear similar in decision space (Figs. 6(b)-(c)), Neural-ES clearly outperforms in objective space, with most EPSL-R1 solutions being dominated by Neural-ES solutions.\n\n------ \n\n> ### **W3**: Results on the trajectory planning problem\n\n**A3:** We sincerely appreciate your valuable feedback regarding both the real-world application section and the presentation of results. Below, we provide both quantitative and additional experimental evidence to support our claims: \n\n| | Neural-ES | EPSL-R1 | EPSL |\n|:---------------------:|:---------:|:-------:|:-------:|\n| #Solutions dominated by Neural-ES | - | 194/300 | 299/300 |\n| HV (higher the better)      | **8.314** | 5.144 | 4.186 |\n\nWe acknowledge that, according to Fig. 7, EPSL-R1 and Neural-ES possibly excel at different objectives. Nonetheless, because a non-dominated solution is always preferred over a dominated one, and HV usually directly reflects both convergence and diversity, we believe the results objectively demonstrate the efficacy of Neural-ES. \n\n__New experiments on UAV navigation__\nTo further strengthen our claims, we conducted additional experiments on real-world tri-objective unmanned aerial vehicle (UAV) navigation problems with varying difficulty levels, with the following acquired results: \n\n| Difficulty | Neural-ES (HV) | EPSL (HV) | EPSL-R1 (HV) |\n|:----------:|:--------------:|:---------:|:------------:|\n| Easy | **3.397±0.279** | 2.471±0.356 | 1.577±0.149 |\n| Moderate | **3.443±0.452** | 2.399±0.372 | 1.520±0.267 |\n| Hard | **2.442±1.143** | 0.716±0.197 | 0.842±0.139 |\n\nThe performance lead by Neural-ES across all difficulty levels validates its superiority in practical applications.\n\n__Availability of Raw Data__: We sincerely apologize for any confusion regarding the presentation format. Due to the space limit, we included representative visual results in the main text, while attaching in the supplementary document all raw quantitative data for producing Fig.4 (i.e., Table 2), as well as the convergence curves for all instances (i.e., Fig. 15). \n\n------\n\n> ### **Q2**: Details about IGD and HV\n\n**A4:** Given a subset of the ground-truth Pareto front $\\widehat{PF}$, and an approximate Pareto front $PF$ by an algorithm, the IGD value is given by\n\n$IGD(PF,\\widehat{PF})=\\frac{1}{|\\widehat{PF}|}\\sum\\limits_{F^\\star\\in\\widehat{PF}}\\min\\limits_{F\\in PF}\\|F-F^\\star\\|_2$.\n\nIntuitively, IGD measures the distance from each point $F^\\star$ in the ideal $\\widehat{PF}$, to the closest $F$ in $PF$. When imposing ideal properties on $\\widehat{PF}$, the IGD is minimized only when the approximation set $PF$ spreads widely and evenly, or in other words, it converges to a good set. Otherwise, there will always be some $F^\\star\\in\\widehat{PF}$ that remains distant from $PF$, resulting in a large IGD. \n\nThe HV measures the volume of a specific region in the objective space, which is between a solution set and a reference point. It is given by: \n\n$HV_u(S)=\\Lambda(\\{y|\\exists x\\in S: F(x)\\prec y\\prec u\\})$,\n\nwhere $S$ is a solution set, $u\\in\\mathbb{R}^M$ is a reference point in the $F$-space, $\\Lambda$ denotes the Lebesgue measure, and $\\prec$ denotes Pareto dominance. \nThe volume always increases when the solution set approaches more closely to the ground truth PF or spreads more widely across the PF, as they both stretch the region.\n\nBoth HV and IGD measure the convergence and diversity of a solution set to the ground-truth PF. IGD is favoured for its computational efficiency, while HV exhibits preferable properties such as submodularity and monotonicity [1,2]. For years, these two metrics have been widely adopted in the literature. \n\n------\n\n> ### **Limitations**: \n\n**A5:** We appreciate your perspective on the scope of our limitation discussion. Constraints have consistently been a long-standing issue for both multiobjective optimization and evolution strategies [3,4], and we argue that it represents a fundamental technical barrier specific to our Neural-ES approach's current capabilities. Below, please see our detailed justification. \n\n__Manifold discontinuity__: Constrained MOPs often produce fragmented Pareto fronts [5]. This fundamentally breaks our current continuous mapping assumption between preference vectors and solutions.\n\n__Sampling breakdown__: The sampling of a random solution $x=(x_1^T,…,x_p^T )^T\\sim\\mathcal{N}(\\varphi_\\theta(x_{0:i-1}))$ becomes problematic when early segments $x_0$ to $x_{i-1}$ are frequently infeasible, causing cascading failures in generating valid $x_i,...,x_p$. \n\n------\n\n> ### **References** \n\n[1] AP Guerreiro, CM Fonseca, L Paquete. \"The hypervolume indicator: Computational problems and algorithms.\" ACM Computing Surveys (CSUR) 54.6 (2021): 1-42.\n\n[2] JG Falcón-Cardona, CAC Coello. \"Indicator-based multi-objective evolutionary algorithms: A comprehensive survey.\" ACM Computing Surveys (CSUR) 53.2 (2020): 1-35.\n\n[3] D., Paul, and N. Hansen. \"Augmented Lagrangian, penalty techniques and surrogate modeling for constrained optimization with CMA-ES.\" Proceedings of the Genetic and Evolutionary Computation Conference. 2021.\n\n[4] X. Chu,  M. Fei, and W. Gong. \"Competitive multitasking for computational resource allocation in evolutionary constrained multi-objective optimization.\" IEEE Transactions on Evolutionary Computation (2024). \n\n[5] Z. Ma and Y. Wang. Evolutionary Constrained Multiobjective Optimization: Test Suite Construction and Performance Comparisons. IEEE Transactions on Evolutionary Computation, 23(6):972–986, December 2019."}, {"review_id": "WGEnfBLt7k", "reviewer": "Reviewer_W93V", "summary": "This paper proposes Neural-ES, a black-box optimization algorithm for Pareto set learning in high-dimensional multi-objective problems. The method replaces explicit full covariance matrices with a neural network that outputs a low-rank covariance structure in an autoregressive sampling framework. The authors demonstrate the algorithm on the proposed BBO-PSL benchmark suite and a real-world trajectory planning problem, showing competitive or superior performance compared to state-of-the-art MOEAs.", "questions": "1. Could the authors clarify why the results for EPSL-R1 are not reported in Figure 6(d)? \n\n2. Could the authors provide more details on how the BBO-PSL benchmark suite was designed and constructed?", "limitations": "Yes", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths: \n\n1. The paper clearly articulates why existing methods fail to scale in high-dimensional multi-objective settings. The problem is well-motivated and practically relevant. \n\n2. The mathematical derivations, including the gradient estimation and the natural gradient correction, are rigorous and clearly presented. \n\n3. The authors designed the BBO-PSL benchmark, which covers a diverse set of challenging properties (nonseparability, scalability, known Pareto sets). \n\n\nWeaknesses:\n\n 1. The hyperparameters are fixed across experiments, without a sensitivity analysis. \n\n2. The paper does not include sufficient ablation experiments. This makes it harder to attribute performance gains to specific components. \n\n3. No runtime measurements are reported. It is unclear how much additional time the autoregressive sampling introduces compared to standard ES variants. This makes it hard to assess the method’s practical efficiency.\n\n 4. Most algorithmic components (low-rank approximation, natural gradient) are adaptations of prior work. The main novelty is in combining autoregressive sampling with ES.", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely appreciate the reviewer's thorough evaluation and constructive feedback. Below, we provide point-by-point responses to address all raised concerns with a clear distinction between:\n* Existing content (simple past tense)\n* New analyses (present perfect tense)\n\n------\n\n> ### **W1:** Hyperparameter sensitivity analysis\n\n**A1:** As detailed in Section C.4 of the supplementary material, we conducted sensitivity analyses on two critical hyperparameters, which are: \n* $p$, the partition count, and \n* $N$, the number of sampled subproblems. \n\nAs mentioned, the results in Fig. 16 indicated that \n* A larger $p$ usually leads to better performance, because it increases the number of Gaussian components and thus enhances the ability to decouple dimensional dependencies; \n* Either too large or too small of $N$ deteriorates the performance, because the former consumes too much evaluation budget while the latter reduces the precision of the estimated gradient. \n\nThroughout the paper, we set $p=n/16$ as it balances performance and computational efficiency, and we set $N=5n/16$, which yields the best performance in general, according to Figs. 16(c) and (d). \n\n------\n\n> ### **W2:** Ablation study\n\n**A2:** Neural-ES is designed based on a minimalistic principle, containing few modular components. As illustrated in Fig.1, Neural-ES primarily consists of two parts: the set model and the back-end search distribution. Our core contribution is a novel formulation of the search distribution, which endows PSL with the ability to effectively deal with dimensional dependencies. \n\nThe experiments on EPSL and EPSL-R1 in Sec. 5 serve dual purposes: they are both comparisons between Neural-ES and the state-of-the-art methods and ablation studies. More specifically, \n* The comparison versus EPSL-R1, which is synthesized by “setting p = 1 in Neural-ES”, and “represents the simplest way to integrate PSL with an ES” (lines 276~277), validates the proposed search distribution;\n* The comparison versus EPSL, which is “a finite difference approach is employed to train the set model” (line 113), validates the ability to handle dimensional dependencies. \n\n------\n\n> ### **W3:** Runtime benchmarking\n\n**A3:** Acknowledging the necessity for runtime assessment, we have supplemented new experiments that will be included in the revised manuscript. We have measured wall-clock time on the 1024-dimensional F1/F5 instances of our proposed BBO-PSL suite. The evaluation budget is 4000n for F1, and 2000n for F5. The results are reported as follows: \n| Algorithm | F1 | F5 |\n|:------------:|:---------:|:---------:|\n| LMO-CSO | 1019s (16.9min) | 1432s (23.9min) |\n| IM-MOEA/D | 1044s (17.4min) | 911s (15.2min) |\n| EPSL | 486s (7.76min)  | 308s (5.14min) |\n| EPSL-R1 | 731s (12.2min)  | 443s (7.39min) |\n| Neural-ES | 1042s (17.4min) | 567s (9.45min) |\n\nComparing the PSL algorithms with the traditional MOEAs, we can observe that the former run faster, or at least as fast, as the latter. This is because most MOEAs rely on additional techniques to maintain population diversity and to deal with slow convergence. Most, if not all, of these techniques come with high overhead and scale poorly with increasing population sizes. In contrast, the PSL paradigm does not rely on any of these techniques, and it turns out that training the set model is less costly than usually believed when compared to MOEAs. \n\nComparing Neural-ES to EPSL and EPSL-R1, the former runs slightly slower than the latter two, but remains competitive overall. As mentioned in Sec.B7 of the supplementary document, Neural-ES yields $\\mathcal{O}(n)$ time complexity when ignoring the network inference cost and $\\mathcal{O}(n^2)$ when including it. To compare, EPSL and EPSL-R1 exhibit the same complexities, only that they have lower overhead due to simpler search distributions. \n\n------\n\n> ### **Q1:** Display issue\n\n**A4:** In Fig. 6(d), solutions generated by EPSL-R1 were unreported because they were dominated by solutions of Neural-ES and EPSL. As noted in the caption of Fig.6, we only displayed solutions that were not dominated by others for clarity. Similar situations occurred in Figs. 9(e) to (h) of the supplementary document, where, at times, the results for EPSL-R1 were not displayed, while at other times, it could be the other algorithms that were missing. \n\n------\n\n> ### **Q2:** Details of the BBO-PSL suite\n\n**A5:** The details of the BBO-PSL suite were presented in Sec. C.1 of the supplementary document. In brief, the instances in the suite were designed based on Theorem 1 in the supplementary document. Three functions $f^{pf}$, $f^{ps}$, and $f^{dis}$ independently control the shapes of the PF and PS, as well as the optimization difficulty. Notably, the distance function $f^{dis}$ injects complex dimensional dependencies into the problems by utilizing the popular rotated ellipsoidal function (i.e., Eq. (31)). We depict the ground-truth PSs in the ($x_1-x_2-x_3$) subspace in Fig. 8, which shows their diverse and highly non-linear structures."}, {"review_id": "DOKQV1I0A7", "reviewer": "Reviewer_evfh", "summary": "The author proposes a new method called Neural-ES, which integrates evolutionary strategies into the paradigm of Pareto set learning. The key idea is to use neural networks to model and capture the complex dimensional dependencies in high-dimensional design variables, rather than relying on the traditional evolution strategy covariance matrix.", "questions": "See Weaknesses part.", "limitations": "Yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "Strengths:\n1. This paper proposed a novel method for solving the Black-box Pareto Set Learning problem.\n\n2. The experimental setting is perfect, and a large number of numerical test results prove the effectiveness of the method.\n\nWeaknesses:\n1. I checked the provided three refs on Line 230. These papers do not seem to provide a clear convergence rate of the ES method. The author should provide more discussion on the theoretical convergence properties of this method. Can it directly link to the convergence analysis of the standard ES method? Can the authors add a discussion with \"Black-box Optimizer with Implicit Natural Gradient\" about the connection between methods and convergence?\n\n2. I still have concerns about the Pareto properties of the solved MOO. It seems the original problem is solving $F(x)$. But the proposed method is finding the gradient of $J(\\theta)$, is the solution of optimizing $J(\\theta)$ the solution of minimizing $F(x)$. In the convex case, I think the ES method works, but in the non-convex case, the Pareto-stationay point of the distribution approximated objective and original objective seems not equal. (See \"Adaptive Stochastic Gradient Algorithm for Black-box Multi-Objective Learning\")\n\n3. Is there any practical application scenario here, about Pareto Set Learning in the black box scenario? We have lots of Pareto optimization methods, such as PaMaL, LORPMAN, COSMOS, PHN, why do we need it (for the black-box case). In many numerical cases, I think the advantage of black-box optimization is it can solve non-smooth or non-continuous, or complex functions, but ES-based method seems do not have a theoretical guarantee on those problems.", "was_revised": true, "final_justification": "My concerns have been addressed, so I maintain my positive score.", "rebuttal": "Thank you for your insightful and constructive comments and suggestions. We appreciate that you recognized our method as novel and commended our experimental setting. Below, please find our detailed responses. \n\n------\n\n> ### **W1** Theoretical convergence rate\n\n**A1**: We appreciate the insightful and constructive comments and break down our response into the following aspects: \n\n__Analyses on the convergence rate of ES remain challenging__. Recently have seen a significant development towards the convergence rate of ES methods, where some empirical and theoretical evidence suggests that ES converges linearly on some specific functions. For example, both Glasmachers [1] and Akimoto et al. [2] utilized the drift analysis to prove the convergence, where an overall Lyapunov function of the state of the algorithm (distribution mean and step-size) is used to prove upper and lower bounds of the expected first hitting time of an epsilon neighborhood of the optimum. However, most of the existing works heavily rely on fierce assumptions, which can be unfounded in practice. For example, the analyses of Glasmachers and Akimoto are limited to the simple (1+1)-ES, and the scale of the covariance matrix in Akimoto’s method needs to be strictly bounded. \n\n__Convergence analysis for multiobjective optimization algorithms is limited__. In multiobjective optimization, seeking a solution set rather than a single Pareto optimal solution makes the analysis even more challenging. This is because the most commonly used Pareto dominance is a partial order rather than a total order, making it difficult to define the optimality for a set-based optimization algorithm. Although analyzing the convergence to a single Pareto stationary point remains possible, which has been elegantly pursued in works like [3], there has been limited research, to the best of our knowledge, investigating the theoretical convergence under Pareto dominance. \n\n__Discussions with INGO__: In [4], Lyu and Tsang proposed a black-box optimization algorithm called Implicit Natural Gradient Optimization (INGO). To avoid explicitly inverting the Fisher, they formulated a surrogate objective where the KL-divergence is attached as a penalty term and can be explicitly derived (Theorem 3). The authors showed that the gradient with respect to the expectation of the sufficient statistics amounts to the natural gradient of the natural parameter of an exponential-family distribution. Furthermore, the authors proved that INGO achieves sublinear convergence when the surrogate objective is $\\gamma$-strong convex (Theorem 4). For Lipschitz continuous, convex, and black-box functions, the authors proved that INGO also achieves sublinear convergence, the overhead of which can be further improved when applying orthogonal sampling (Theorem 5). We compare Neural-ES with INGO as follows: \n* Both Neural-ES and INGO avoid explicitly inverting the Fisher, which substantially reduces the variance of the estimated gradient and increases the computational efficiency. \n* Both the empirical evidence in Fig. 15 of the supplementary document and the theoretical results of INGO suggest that both algorithms achieve sublinear convergence. \n* It may be feasible to replace Eq. (10) with the implicit natural gradient formulation of INGO, but the subsequent deviations need to be changed accordingly. \n* However, as mentioned above, the primary difficulty of convergence analysis with the PSL paradigm arises from the lack of a definition of optimality based on a total order. As a result, the convergence analysis used in INGO may not be easily applicable to Neural-ES. \n\nTo summarize, recent findings show that some ES algorithms achieve linear global convergence, but strong assumptions are required; meanwhile, the difficulty of obtaining the convergence rate for set-based multiobjective optimization algorithms (and PSL as well) mostly arises from the lack of a clear definition of total-order optimality. These two main reasons obstruct the analysis for the proposed Neural-ES. Nonetheless, we believe your insightful suggestion clearly helps improve the quality of our work, and will add the above discussions to the revised manuscript. \n\n------\n\n### **W2:** Convergence of the surrogate $\\mathcal{J}(\\theta)$\n\n**A2:** In [3], Ye et al. proposed the Adaptive Stochastic Multi-objective Gradient (ASMG) algorithm for black-box MOPs, extending the INGO algorithm with promising results. A key finding is that, in black-box scenarios, ASMG may only converge to an $\\epsilon$-accurate Pareto stationary point rather than an exact one (Propositions 4.5, 4.6). From our understanding, this is because the Gaussian smoothing relaxes the original KKT conditions in MOPs, which means the gradient norm of the original $F$ can still be slightly greater than zero, but is still strictly less than the upper bound set by the covariance norm, when the surrogate $\\mathcal{J}$ achieves Pareto stationary. We fully agree with this observation as it may be a common issue for most black-box MOO algorithms, including Neural-ES. We appreciate you bringing it to our attention, and will clarify this point in detail (e.g., in both Sec. 3.1 and the limitation) in the revised manuscript.\n\n------\n\n> ### **W3**: Practical applications\n\n**A3**: Fully acknowledging the necessity for more practical applications, we have supplemented our work with new experiments on a tri-objective unmanned aerial vehicle (MO-UAV) navigation problem. The problem is included in the Meta-BBO-v2 library [5,6]. We have examined the PSL algorithms on three problem instances with varying difficulty levels. The primary results are as follows: \n\n| Instance | Dimensions | Neural-ES (HV) | EPSL (HV) | EPSL-R1 (HV) |\n|:---------:|:----------:|:--------------:|:------------:|:-------------:|\n| Easy | 30 | **3.397±0.279**  | 2.471±0.356 | 1.577±0.149 |\n| Moderate | 60 | **3.443±0.452**  | 2.399±0.372 | 1.520±0.267 |\n| Hard | 120 | **2.442±1.143**  | 0.716±0.197 | 0.842±0.139 |\n\nThe results highlight the efficiency of Neural-ES. On one hand, PSL is capable of most categories of black-box MOPs as most other algorithms (e.g., MOEAS) do; on the other hand, its merit lies in the ability to \"provide a tailored optimal solution for any valid decision-maker preference trade-off without requiring re-optimization from scratch\", as we stated in lines 31-32 of the main text. \n\nBesides, we totally agree with you that the \"advantage of black-box optimization is it can solve non-smooth or non-continuous, or complex functions\". In fact, some recent studies [7,8] have proven that ES converges globally and linearly, on more than convex functions, such as quasiconvex functions and discontinuous functions, as well as those with non-convex sublevels. Therefore, combining ES with PSL could be one of the many promising directions for tackling black-box MOPs. \n\n------\n\n> ### **References**\n\n[1] T. Glasmachers. \"Global convergence of the (1+1) evolution strategy to a critical point.\" Evolutionary computation 28.1 (2020): 27-53. \n\n[2] Y. Akimoto, et al. \"Global linear convergence of evolution strategies on more than smooth strongly convex functions.\" SIAM Journal on Optimization 32.2 (2022): 1402-1429.\n\n[3] F. Ye, et al. \"Adaptive stochastic gradient algorithm for black-box multi-objective learning.\" ICLR. 2024.\n\n[4] Y. Lyu, and I.W. Tsang. \"Black-box optimizer with stochastic implicit natural gradient.\" ECML-PKDD, 2021.\n\n[5] Z. Ma, et al. \"Metabox: A benchmark platform for meta-black-box optimization with reinforcement learning.\" NeurIPS 36 (2023): 10775-10795.\n\n[6] Z. Ma, et al. \"MetaBox-v2: A Unified Benchmark Platform for Meta-Black-Box Optimization.\" arXiv preprint arXiv:2505.17745 (2025).\n\n[7] A. Auger, and N. Hansen. \"Linear convergence of comparison-based step-size adaptive randomized search via stability of Markov chains.\" SIAM Journal on Optimization 26.3 (2016): 1589-1624. \n\n[8] C. Toure, A. Auger, and N. Hansen. \"Global linear convergence of evolution strategies with recombination on scaling-invariant functions.\" Journal of Global Optimization 86.1 (2023): 163-203."}]

Author rebuttal:

Meta-review: This paper proposes a neural evolution strategy for high-dimensional black-box multi-objective optimization problems via Pareto set learning (PSL). The Gaussian distributions modeled by a neural network are used for solution generation. The neural network is updated using the estimated gradient by leveraging the natural gradient of the Gaussian distribution. The proposed Neural-ES is evaluated using a benchmark suite for black-box PSL. The experimental result shows the efficiency of Neural-ES, demonstrating its capability to effectively learn the Pareto front in several multi-objective scenarios.

This paper integrates evolution strategies with neural network-based PSL, which is an interesting direction. The experimental results verify the effectiveness of the proposed method.

The authors provided extra experiments in the rebuttal and adequately answered the reviewers' concerns. Based on the reviewers' comments and discussion, the contribution of this paper meets the standards of acceptance.

I would encourage the authors to include additional explanations and experiments discussed in the rebuttal phase in the revised paper. I also think that it is important to justify the reformulation of the original problem using Eq. (4) as pointed out by the reviewer.

In addition, as a minor comment from AC, I could not find the setting of the learning rate $\eta$ for the proposed method. It would be better to clarify it and analyze the sensitivity of the learning rate setting. Also, in the rebuttal, the authors mentioned that the proposed method used the Adam optimizer, but I could not find its description in the main paper. It would be better to clarify it in the revised paper.
Decision: Accept (poster)