LibMOON: A Gradient-based MultiObjective OptimizatioN Library in PyTorch
Xiaoyuan Zhang Liang Zhao Yingying Yu Xi Lin Zhenkun Wang Han Zhao Qingfu Zhang

Continual Diffuser (CoD): Mastering Continual Offline Reinforcement Learning with Experience Rehearsal
      Dacheng Tao^8
    3 September 2024
======================================================================================================

§ ABSTRACT 
    Multiobjective optimization problems (MOPs) are prevalent in machine learning, with applications in multi-task learning, learning under fairness or robustness constraints, etc. Instead of reducing multiple objective functions into a scalar objective, MOPs aim to optimize for the so-called Pareto optimality or Pareto set learning, which involves optimizing more than one objective function simultaneously, over models with millions of parameters. Existing benchmark libraries for MOPs mainly focus on evolutionary algorithms, most of which are zeroth-order methods that do not effectively utilize higher-order information from objectives and cannot scale to large-scale models with millions of parameters. In light of the above gap, this paper introduces , the first multiobjective optimization library that supports state-of-the-art gradient-based methods, provides a fair benchmark, and is open-sourced for the community.






§ INTRODUCTION

MultiObjective Optimization problems (MOPs) are ubiquitous in machine learning. For instance, trustworthy machine learning (TML) includes fairness problems balancing the fairness level and the accuracy level <cit.>; 
in robotics, it is necessary to balance several objectives (e.g., forward speed and energy consumption)  <cit.>; similarly, recommendation systems face potentially conflicting objectives, such as novelty, accuracy, and diversity <cit.>. 
For all the applications above, the underlying optimization problem involves an MOP with m objectives and can be (informally) defined as:

    min_θ(θ) = ( L_1(θ), …, L_m(θ) ),

where L_1(θ), …, L_m(θ) denote m (potentially) conflicting objectives and we denote the size of the model parameter as N  |θ|. Note that as informally defined above, <Ref> is a vector optimization problem that does not necessarily admit a total ordering. Nevertheless, we use the formulation above to emphasize that the underlying problem involves m objectives and we will formally define the goal as follows.

For non-trival multiobjective problems, no single solution can attain the minimum of all objectives simultaneously. To compare two solutions for an MOP, we introduce the concepts of dominance and Pareto optimality <cit.>. 
Dominance occurs when solution θ^(a) satisfies L_i(θ^(a)) ≤ L_i(θ^(b)) for all 1 ≤ i ≤ m, with at least one strict inequality L_j(θ^(a)) < L_j(θ^(b)).
A solution is Pareto optimal if no other solution in the feasible region dominates it. The set of all Pareto optimal solutions is the Pareto set (PS), and its image is the Pareto front (PF).

Over the last few decades, multiobjective evolutionary algorithms (MOEAs) emerged as a widely used methodology for addressing MOPs. 
Due to their population-level designs, MOEAs can identify a diverse set of Pareto optimal solutions that approximate and represent the entire PF of a given MOP.
Several libraries have emerged as standard platforms for fair comparison of MOEAs, including PlatEMO  <cit.> (in Matlab supporting over 160 MOEAs), Pagmo  <cit.> (in C++), and Pymoo  <cit.> (in Python). 

Compared to MOEAs, gradient-based multiobjective optimization (MOO) methods are particularly suitable designed for large-scale machine learning tasks involving thousands to millions of neural network parameters. While these methods can find Pareto stationary solutions—solutions that cannot be locally improved in any objective—they are also effective approximators of Pareto solutions for nonconvex models, including deep neural networks.

In this paper, we roughly divide recent gradient-based MOO methods into two main groups.
The first one focuses on finding a finite set of Pareto solutions to approximately represent the entire Pareto front, exemplified by the
Upper Bounded Multiple Gradient Descent Algorithm (MGDA-UB) <cit.> and Exact Pareto Optimization (EPO) <cit.>. The second, known as 

Pareto set learning (PSL) <cit.>, develops a single neural model parameterized by ϕ to represent the entire Pareto set. Since the parameter space of ϕ is large, thus PSL is usually optimized by gradient-based methods.   

With the growing use of gradient-based MOO methods in neural networks <cit.>, there is a pressing need for a standard library to benchmark comparisons and developments. To fulfill this demand, we introduce , a gradient-based Multi-Objective Optimization Library supporting over twenty state-of-the-art methods.

We summarize our contributions as follows:

    
  * We propose the first modern gradient-based MOO library, called .  is implemented under the PyTorch framework <cit.> and carefully designed to support GPU acceleration.  supports synthetic problems and a wide range of real-world multiobjective machine learning tasks (e.g., fair classification and multiobjective classification problems).
    
  *  supports over twenty state-of-the-art (SOTA) gradient-based MOO methods for constructing Pareto optimal solution sets, including those using finite solutions to approximate the whole PS/PF <cit.>, PSL methods <cit.> aimed at approximating the entire Pareto set with a single neural network, and MOBO methods that designed to avoid frequent function evaluations. 
    
  *  has already been open-sourced at <https://github.com/xzhang2523/libmoon> with document at <https://libmoondocs.readthedocs.io/en/latest/>, offering a reliable platform for comparing SOTA gradient-based MOO approaches and providing standard multiobjective testing benchmarks. Beyond examining its source code,  can be installed via  as an off-the-shelf gradient-based multiobjective tool.




  
Notation. In this paper, bold letters represent vectors (e.g.,  for preference vectors, (·) for vector objective functions), while non-bold letters represent scalars or network parameters (e.g., θ). The number of objectives is denoted by m, and K represents the size of finite solution set. The preference vector  lies in the (m-1)-dimensional simplex, satisfying ∑_i=1^m λ_i=1 and λ_i ≥ 0. The decision network parameter θ has a size of n. Refer to <Ref> for the full notation table.




§ RELATED WORKS

In this section, we review some methods for gradient-based multiobjective optimization in <Ref> and existing multiobjective optimization libraries in <Ref>.



 §.§ Gradient-based multiobjective optimization
 
Gradient-based MOO has a long research history. For instance, the famous convex optimization book <cit.>Chap 4 discusses using linear scalarization to convert MOO into a single-objective optimization (SOO) problem. Historically, gradient-based methods were not mainstream for multi-objective optimization, with evolutionary algorithms (MOEAs) being more prevalent due to their ability to maintain a population, aligning well with the Pareto set. However, gradient-based MOO has gained popularity recently with advancements in machine learning. A seminal work is MGDA-UB <cit.>, which introduced MOO concepts into deep learning by framing MTL as MOO. Since then, numerous approaches like EPO <cit.>, Pareto MultiTask learning (PMTL) <cit.>, and MOO with Stein Variational Gradient Descent (SVGD) <cit.>, including methods for learning the entire Pareto set <cit.>, have emerged. This paper aims to integrate these developments into a unified framework.




 §.§ Multiobjective optimization libraries
 
There a number of multiobjective library before our work. The difference between  and previous libraries is that is designed to support gradient-based MOO methods. Comparison can be found in <Ref>. 

LibMTL. <cit.> LibMTL is a Python-based library for multitask learning. The key difference between and LibMTL is that LibMTL focuses on finding a single network to benefit all tasks, such as optimizing direction or network architecture. In contrast,  addresses inherent trade-offs, where improving one objective inevitably worsens others, and explores the distribution of Pareto solutions.

jMetal. <cit.> jMetal is a comprehensive Java framework designed for multi-objective optimization, leveraging the power of metaheuristics. It offers a highly flexible and extensible platform, making it accessible for users across various domains. With its user-friendly interface and robust architecture, jMetal has become a popular choice for researchers and practitioners in diverse application areas.


Pymoo. <cit.>
Pymoo is a Python-based framework for multiobjective optimization problems with many local optima. It supports mainstream MOEA methods such as NSGA-II <cit.>, NSGA-III <cit.>, MOEA/D <cit.>, and SMS-EMOA <cit.>. Pymoo allows flexible algorithm customization with user-defined operators and includes tools for exploratory analysis and data visualization. It also offers multi-criteria decision-making strategies to help users select optimal solutions from the Pareto set.

PlatEMO <cit.>.
PlatEMO is a MATLAB-based multiobjective optimization tool supporting over 160 evolutionary algorithms and a comprehensive set of test problems. It supports diverse problem types, including sparse, high-cost, large-scale, and multimodal. PlatEMO includes performance metrics for evaluating optimization efficacy and offers tools for visualization and interaction during the optimization process.

Pagmo <cit.>. 
Pagmo, a C++ library, executes massively parallel multiobjective global optimization using efficient evolutionary algorithms and gradient-based techniques like simplex, SQP, and interior-point methods. It supports concurrent algorithm execution and inter-algorithm collaboration through asynchronous generalized island models, addressing a wide range of optimization challenges, including constrained, unconstrained, single-objective, multiobjective, continuous, integer, stochastic, and deterministic problems.

EvoTorch <cit.> and EvoX <cit.>.
EvoTorch accelerates evolutionary algorithms in PyTorch, while EvoX scales them to large, complex problems with GPU-accelerated parallel execution for single and multiobjective tasks, including synthetic and reinforcement learning. These libraries highlight a trend toward using GPU parallelization to manage larger problem scales and decision variables, which is also of the interest of this paper. This work focuses on creating a platform for efficient gradient-based multiobjective methods to solve large-scale machine learning problems in PyTorch.






§ : A GRADIENT-BASED MOO LIBRARY IN PYTORCH
 
To introduce  , we first introduce its framework in <Ref>, and briefly introducing its supporting problems and metrics. Then introduce supported solvers in <Ref>.





 §.§ Framework
 
<Ref> demonstrates the framework of the proposed  library.  supports three categories of solvers: MOO solvers aiming to find a finite set of Pareto solutions satisfying certain requirements, PSL solvers aiming to learn whole PS with a single model, and MOBO solvers aiming to solve expensive problems. Each solver category is designed highly modulized and new solvers are easy to plugin  by rewriting only a small portion of code, e.g., the way of gradient manipulation [An example of adding a new method is given in <https://libmoondocs.readthedocs.io/en/latest/develop/add_method.html>.]. MOO and PSL solvers support all synthetic, MTL, and realworld (RE) problems, while MOBO solvers support synthetic and RE problems. 








  
Supported problems.
 currently supports three categories of methods, synthetic problems, MTL problems, and RE problems.



  
Supported metrics.
 supports a wide range of metrics namely, (1) hypervolume (HV), (2) inverted general distance (IGD), (3) fill distance (FD), minimal distance (l_min), (4) smooth minimal distance (sl_min), (5) Spacing, (6) Span, (7) penalty-based intersection (PBI), (8) inner product (IP), (9) cross angle (ϑ). Full descriptions are provided in <Ref>. 




 §.§ MOO solvers
 





  
Aggregation function.
To find Pareto solutions, a straightforward way is to use an aggregation function g_λ^agg : ℝ^m →ℝ to transform a multi-objective problem (MOP) into a single-objective optimization problem (SOP). 
This involves minimizing:

    min_θ g_λ^agg((θ)),

where (θ) represents the objective vector and g_λ^agg is an aggregation function determined by preference λ. If g_λ^agg is decreasing with respect to (θ) (i.e., g_λ^agg((θ^(a))) < g_λ^agg((θ^(b))) when _i(θ^(a)) ≤_i(θ^(b)) for all i with at least one strict inequality), the optimal solution of <Ref> corresponds to a Pareto optimal solution of the original MOP (<Ref>) <cit.>.

The linear scalarization (LS) function, g_λ^agg((θ)) = ∑_i=1^m λ_i L_i(θ), is the most commonly used aggregation function. Other notable functions include the (Smooth) Tchebycheff (Tche) <cit.>, modified Tchebycheff (mTche) <cit.>, Penalty-Based Intersection (PBI) <cit.>, and COSMOS <cit.>. For detailed formulations, see <Ref>. Aggregation-based methods can be efficiently optimized in  using single objective backpropagation. Aggregation-based method update the parameter θ iteratively as follows:

    θθ - η = θ - η∂ g_λ^agg((θ))/∂θ,

where  is the updating direction.  



  
Generalized aggregation function. 
In addition to optimizing the single-objective aggregation function, multiple gradient manipulation methods adjust gradients and compute update directions to ensure solutions satidfy specific properties. 
These methods can be summarized as follows: the update direction  is a dynamically weighted sum of gradient vectors <cit.>, expressed as = ∑_i=1^m α̃_i ·∇ L_i (θ).

Equivalently, updating the direction involves optimizing a generalized aggregation function.

    g̃_λ((θ)) = α̃_i · L_i(θ).




    Multiple gradient manipulation methods.
 


As listed in <Ref>, notable methods include EPO <cit.> for finding `exact Pareto solutions' (intersection points of Pareto front and preference vectors), HVGrad <cit.> for maximizing hypervolume, MGDA-UB <cit.> for finding directions beneficial to all objectives, MOO-SVGD <cit.> for diverse solutions, PMGDA <cit.> for Pareto solutions satisfying specific user requirement, and PMTL <cit.> for sector-constrained Pareto solutions. 




In  , multiple gradient manipulation methods are implemented by first calculating dynamic weight factors α̃_i and then performing backpropagation on the generalized function g̃_λ((θ)) to update the neural network parameters. Those weight factors α̃_i are typically computed by solving a linear programming (LP) problem (e.g., <cit.>[Eq. 24], <cit.>[Eq. 16]), a quadratic programming (QP) problem (e.g., <cit.>[Eq. 14], <cit.>[Eq. 3]), or through other methods like hypervolume gradient estimation <cit.>.

Some methods use user-specific preference vectors as input, making the aggregation function dependent on preference , while others do not (preference-free optimizers). A summary of whether an optimizer is preference-based or preference-free is provided in <Ref>.



    Zero-order optimization.
 
 not only supports first-order optimization when the gradients of those objective functions are known, but also supports zero-order optimization methods with estimated gradients ∇̂ L_i(θ) using methods such evolutionary strategy (ES) <cit.>. 




 §.§ Pareto set learning solvers





In addition to the existing methods for finite Pareto solutions,  supports Pareto Set Learning (PSL), which trains a model with parameter ϕ to generate an infinite number of approximate Pareto solutions.
We denote the model output as θ and represent the model as θ_ϕ(·) : Δ_m →ℝ^n, where Δ_m is the m-dim preference simplex and ℝ^n is the entire decision space of θ.

PSL can be deemed as learning a model that find solutions under infinitely many preference well. 
The loss function of PSL is formulated as follows:  

    min_ϕℓ_psl = _∼Unif(Δ_m)g̃_(( θ_ϕ() )),     θ_ϕ(·) : Δ_m ^n.

The function g̃_(·): ℝ^m →ℝ is a generalized aggregation function (as introduced in <Ref>), which converts an m-dimensional objective vector into a scalar based on specific preferences. Since <Ref> involves preference expectation, all preference-based finite Pareto solvers (<Ref>) can serve as basic Pareto solution solvers, while preference-free solvers are unsuitable for PSL.


For PSL, we explain the function (θ_ϕ()) with more details as follows: In synthetic problems, the decision variable is directly the model output θ_ϕ, and θ_ϕ is evaluated by the loss function. In multi-task learning (MTL) problems, θ_ϕ() represents the parameters of a target network, and (θ_ϕ()) is evaluated using the dataset (x_i, y_i) with this target network. See <Ref> for a visualization of these architectures.

The gradient of ℓ_psl can be estimated by the chain rule:


    ∂ℓ_psl/∂ϕ_1 × D = _∼Unif(Δ_m)g̃__A: (1 × m)·θ_B: (m × n)θϕ_C: (n × D).



The expectation is empirically estimated using a batch of K preferences. The gradient vector A=(α̃_1, …, α̃_m) is computed by the finite Pareto solution solver. The gradient matrix C can be calculated either via backpropagation (when gradients are easily obtained) or using a zero-order method. In <Ref>, C is always estimated through backpropagation, since θ is a continuous function of ϕ. Gradient calculation by existing PSL methods is summarized in <Ref>. After calculating the gradient of ℓ_psl, the parameter ϕ is iteratively updated as ϕϕ - η∂ℓ_psl/∂ϕ.



 §.§ Multiobjective Bayesian Optimization solvers


When the evaluation of objective functions is costly, multiobjective Bayesian optimization (MOBO) is often the preferred methodology for tackling such challenges. While there are several existing libraries, such as Botorch <cit.> and HEBO <cit.>, that facilite MOBO, they largely overlook algorithms that leverage decomposition techniques like PSL. To bridge this gap,  also includs three decomposition-based MOBO algorithms, including PSL-MOBO <cit.>, PSL-DirHVEI <cit.>, and DirHV-EGO <cit.>. In each iteration, these methods build Gaussian process (GP) models for each objectives and generate a batch of query points for true function evaluations.

PSL-MOBO is an extension of PSL method for expensive MOPs. It optimizes the preference-conditional acquisition function α(θ|) over an infinite number of preference vectors to generate a set of promising candidates: 




    min_ϕℓ_psl = _∼Unif(Δ_m) [α(θ_ϕ()|)],     θ_ϕ() : Δ_m ^N.

 

Currently, our library supports two representative preference-conditional acquisition functions: the Tchebycheff scalarization of lower confidence bound (Tch-LCB) <cit.> and the expected direction-based hypervolume improvement (DirHV-EI) <cit.>. We note that DirHV-EI can be regarded as an unbiased estimation of a weighted expected hypervolume improvement. Our library also supports DirHV-EGO, which employs a finite set of predetermined reference vectors as in <cit.>.



§ SNIPPETS OF  EXAMPLES

For the reader's convenience, we demonstrate some snippets to show how  solves synthetic problems efficiently by just a few lines of codes. For example codes solving multitask learning problems, please refer to <Ref>.  can be installed by .

[style=pythonstyle, caption=Finding a size-K (K=5) Pareto solutions with four lines of code in .] 
from libmoon.solver.gradient.methods import EPOSolver
from libmoon.util import synthetic_init, get_problem, uniform_pref

problem = get_problem(problem_name='ZDT1')
prefs = uniform_pref(n_prob=5, n_obj = problem.n_obj, clip_eps=1e-2)
solver = EPOSolver(problem, step_size=1e-2, n_iter=1000, tol=1e-2)
res = solver.solve(x=synthetic_init(problem, prefs), prefs=prefs)


[style=pythonstyle, caption=Pareto set learning in a problem with three lines of solving problem and two lines of evaluating the results in .] 
from libmoon.solver.psl.core_psl import AggPSLSolver
from libmoon.util import synthetic_init, get_problem, uniform_pref

problem = get_problem(problem_name='ZDT1')
# agg list ['ls', 'tche', 'mtche', 'pbi', ... ]
solver = AggPSLSolver(problem, agg='ls')
model = solver.solve()
prefs = uniform_pref(n_prob=100, n_obj=problem.n_obj, clip_eps=1e-2)
eval_y = problem.evaluate(model(Tensor(prefs).cuda()))


Code in Listing 1 demonstrates solving a synthetic problem efficiently with four lines of  code: 1. Define the problem. 2. Specify the preference vectors. 3. Specify the solver. 4. Obtain the results. The result `res' contains five elements: `x_opt' (optimal Pareto solution), `y_opt' (Pareto objective), `x_history' (decision history), `y_history' (objective history), and `hv_history' (hypervolume history). Codes in Listing 2 show that a Pareto set model can be trained by three lines of codes by assigning a user-specific aggregation function and produce predicted Pareto solutions by two lines of code.



§ EMPIRICAL STUDIES
 
In this section, we study empirical results of . 
Codes are basically run on a person computer with an i7-10700 CPU, and a 3080 GPU with 10G GPU memory.
And we also conduct some analysis regarding to GPU ranging from 3080, 4060, and 4090.
The empirical results contain five parts, 



 §.§ MOO solvers for synthetic problems

We illustrate the performance of various finite Pareto solution solvers using the VLMOP2 problem, formulated as follows:


    min_θ f_1(θ)    = 1 - exp(-∑_i=1^n(θ_i - 1/√(n))^2), 
    min_θ f_2(θ)    = 1 - exp(-∑_i=1^n(θ_i + 1/√(n))^2).

This problem has been widely studied in the literature <cit.> due to its PF has a non-convex shape [In MOO, a PF is convex or non-convex based on whether the objective space is convex or non-convex. The PF represents the non-dominated boundary of the objective space.].  
The visualization results are shown in <Ref> and the numerical results are shown in <Ref>.





We present the following key findings:

    
  * Agg-COSMOS produces approximate `exact' Pareto solutions (the corresponding Pareto objectives aligned with preference vectors). This is due to a cosine similarity term in Agg-COSMOS encouraging the Pareto objectives be close to preference vectors.
    
  * Agg-LS can only find two endpoints on this problem since PF of VLMOP2 is of non-convex shape. Different preference vectors correspond to duplicated Pareto solutions.
    
  * Agg-PBI generates 'exact' Pareto solutions when the coefficient of d_2 exceeds a specific value, as guaranteed by PBI <cit.>. However, this parameter is challenging to tune and is influenced by the Pareto front's curvature. Additionally, PBI can transform a convex multi-objective optimization problem into a non-convex one, making Agg-PBI less recommended.
    
  * Agg-Tche generates diverse solutions and produces exact Pareto solutions corresponding to the inverse of the preference vector. Both Agg-Tche and Agg-SmoothTche retain the convexity of objectives, meaning their aggregation functions remain convex when all objectives are convex. However, the solution distribution of Agg-SmoothTche cannot be precisely determined due to the involvement of a temperature coefficient h.
    
  * HVGrad updates the decision variable using the hypervolume gradient, resulting in the largest hypervolume. PMTL, a two-stage method, can struggle with tuning warm-up iterations, leading to poor initialization and less diverse solutions. MOO-SVGD's update direction has two conflicting goals: promoting diversity and ensuring convergence. This conflict makes MOO-SVGD difficult to converge, often requiring 10 times more iterations than other methods, resulting in poor convergence behavior.
    
  * Among the methods, MGDA-UB, Random, Agg-PBI, and MOO-SVGD exhibit relatively large deviations. MGDA-UB and Random generate arbitrary Pareto solutions due to their computation nature. Agg-PBI results in non-convex aggregation functions, leading to variable solutions based on different random initializations. MOO-SVGD has high standard deviation because its optimization process is unstable, and its balance between convergence and diversity is unclear.
   




 §.§ Pareto set learning on synthetic problems

In this section, we present PSL results using the synthetic VLMOP2 problem as an example. 
Visualizations of the predicted Pareto solutions are shown in <Ref>, and numerical results are reported in <Ref>.



From the table and figure, we have the following key findings:

    
  * PSL with the COSMOS aggregation function fails to find all marginal Pareto solutions because COSMOS does not guarantee the discovery of the entire Pareto set/front. PSL with linear scalarization function could not fit the two endpoints of the PF. 
    
  * PSL with the Smooth Tchebycheff function finds diverse Pareto solutions due to the unclear relationship between preference vectors and Pareto objectives. PSL using Tche., EPO, and PMGDA as base solvers can discover the entire PS/PF, as all three methods effectively find exact Pareto solutions. By traversing the entire preference space, the Pareto model fits the entire Pareto set accurately.  






 §.§ MOO solvers for MTL problems

In this section, we investigate the performance of different finite Pareto solution solvers on the Adult dataset, which is a  multitask learning task called fairness classification. 
For this task, the decision variable θ is the parameter of a fully-connected neural network, with |θ|=28033. 
The first objective function is the cross entropy loss, and the second objective is the hyperbolic tangent relaxation of the difference of equality of opportunity (DEO) loss <cit.>[Eq. 6]. 
The visualization on the Adult problem is shown in <Ref> and the numerical result is shown in <Ref>. 

From the table and figure, we have the following key findings:

    
  * Agg-LS has two main drawbacks: it cannot identify the non-convex parts of the Pareto front (as previously mentioned), and the relationship between preference vectors and Pareto objectives is unknown, making it challenging to select appropriate preference vectors. Agg-PBI and Agg-COSMOS only find a small portion of the PF due to it is hard to set the coefficients.     
    
  * Agg-SmoothmTche and Agg-mTche perform well on this task, as they can accurately find Pareto solutions. Once the Pareto front range is known, diverse solutions can be easily found using uniform preference vectors.
    
  * The Random method and MGDA-UB only find a single Pareto solution, since the position of this solution cannot be controlled by these methods. 
    
  * Among the three methods that directly find a set of Pareto solutions (MOO-SVGD, PMTL, and HV-Grad), HV-Grad produces the most diverse solutions with the largest hypervolume. PMTL, being a two-stage method, can encounter issues when solutions fall outside the sector due to stochastic factors. MOO-SVGD optimizes both convergence and diversity but is generally unstable based on our tests.








 §.§ Pareto set learning on MTL

This section presents the Pareto set learning results on the MO-MNIST problem using a hypernet architecture. The model was trained for 20 epochs, optimizing approximately 3.24M parameters, with the first and second objectives being the cross-entropy losses for the top-right and bottom-left images. EPO-based PSL and PMGDA-based PSL are not very suitable for this task since manipulating gradient of 3.2M dim is not efficient, the visualization results are reported in <Ref> and numerical results in <Ref>.

From the table and figure, we have the following key findings:


    
  * The Pareto Front (PF) of this task is nearly convex, indicating that Local Search (LS) can efficiently find the entire Pareto front. Therefore, it is not necessary to employ complex Multi-Objective Optimization (MOO) methods.
    
  * Agg-LS significantly outperforms other methods on this task, as evidenced by the substantially higher Hypervolume (HV). Furthermore, the training losses across other methods show minimal reduction, indicating that for this challenging task, the convergence of other non-linear scalarization methods is extremely slow.
    
  * Agg-SmoothTche and Agg-Tche produce duplicate marginal Pareto solutions because selecting appropriate weight factors is challenging.












 §.§ MOBO for synthetic and real-world problems

In this section, we test three MOBO algorithms in  on three benchmark problems, including ZDT1, RE21, VLMOP1 and VLMOP2. To ensure a fair comparison, we generate 11n-1 initial samples using Latin Hypercube Sampling for each method. The maximum number of function evaluations is set as 200.
Our experimental results, illustrated in Figure <ref>, clearly demonstrate the rapid convergence capabilities of all three methods. DirHV-EGO, PSL-DirHV-EI, and PSL-MOBO not only efficiently navigate the solution space but also quickly reach optimal solutions. This highlights the robustness and effectiveness of our implemented algorithms in handling different types of MOPs.





 §.§ GPU acceleration

We evaluate  performance on Pareto set learning for the MO-MNIST problem across various platforms (CPU, RTX 3080, 4060, 4090). Running times are detailed in <Ref> and visualized in <Ref>. The table and figure show a significant reduction in time (about one-third) when using a personal GPU compared to a CPU. The RTX 4090 further reduces time by approximately 25% compared to the RTX 4060.








§ CONCLUSION, LIMITATIONS, AND FURTHER WORKS



  
Conclusion.
We introduced the first gradient-based multiobjective optimization framework called , which is implemented based on the PyTorch environment for the research community's convenience. 
 supports mainstream gradient-based MOO methods;
the modular design of  further allows the library to address various MOPs, including discrete Pareto solution optimization, Pareto set learning, multiobjective Bayesian optimization, etc., in a plug-and-play manner.
 can thus be leveraged to quickly yet robustly test new MOO ideas.



  
Limitations. The limitations of  include: (1) the rapid development of gradient-based MOO methods makes it impossible to incorporate all relevant techniques, so some effective methods may be missing; (2) the codebase, primarily based on works from Prof. Qingfu Zhang's group, may not easily accommodate new methods.



  
Future Work. Future work includes: (1) maintaining a user and development community to address issues promptly, and (2) adding newly published or well-known methods to the library as quickly as possible.





§ ACKNOWLEDGE

 was made possible with the valuable feedback and comments from a lot of reseachers. We appreciate the help from Hongzong Li, Xuehai Pan, Zhe Zhao, Meitong Liu, Weiduo Liao, Weiyu Chen, Dake Bu, Yilu Liu, Song Lai, Cheng Gong, Prof. Jingda Deng, Prof. Ke Shang, Prof. Genghui Li, Prof. Zhichao Lu, and Prof. Tao Qin. We sincerely thank them for their early feedback and contributions.




§ APPENDIX
 



 §.§ Full name and notation tables
 
This section lists the full names of optimization methods and terms for clarity (<Ref>) and provides notation in <Ref>. 





 §.§ Aggregation functions

Aggregation function convert an MOP into a single-objective optimization problem under a specific preference vector . Some popular aggregation functions are:




    
  * COSMOS:
    
    g^cosmos(θ) = ^⊤(θ) - μ^⊤(θ)/(θ),

    where μ is a positive weight factor to align the objective vector (θ) with the preference vector .
    
    
  * Linear scalarization (LS):
    
    g^LS(θ) = ∑_i=1^m λ_i L_i(θ).


    
  * Tchebycheff (Tche):
    
    g^Tche(θ) = max_1 ≤ i ≤ m{λ_i (L_i(θ) - z_i) },

    where  is a reference point, usually set as the nadir point the minimal value in each objectives. 

    
  * Modified Tchebycheff (mTche):
    
    g^mTche(θ) = max_1 ≤ i ≤ m{L_i(θ) - z_i/λ_i},

    where z is a reference point, the same as the point used in the Tchebycheff scalarization function. 
    
  * Smooth Tchebycheff (STche):
    
    g^STche(θ) = 1/hlog∑_1 ≤ i ≤ mexp(h ·λ_i (L_i(θ) - z_i)).

    The Smooth Tchebycheff function uses a relaxed Smooth max operator. The advantage of this approach is that g^STche(θ) becomes a Smooth function if each objective function L_i(θ) is Smooth, unlike the non-Smooth g^Tche(θ). Smooth functions generally have faster convergence compared to non-Smooth ones. Similarly, we can define the Smooth modified Tchebycheff function. 

    
  * Penalty-Based Intersection (PBI):
    
    g^PBI(θ) = 
            1/·∑_i=1^m λ_i L_i(θ)_d_1
             + μ(θ) - d_1/·_d_2,

    where μ is a positive weight factor that encourage a objective to align with a preference vector .

    
  * p-norm: 
    
    g^pnorm(θ) = (θ) - _p.

    
  * Augmented Achievement Scalarization Function (AASF):
    
    g^aasf(θ) = g^mtche(θ) + ρ g^LS(θ),

    where ρ is small positive coefficient, usually set as 0.1. 

    The contour curve for this functions for a bi-objective case can be found in <https://libmoondocs.readthedocs.io/en/latest/gallery/aggfuns.html>. 








 §.§ Metrics
 
Metrics can be categorized into two groups. The first group evaluates the quality of a set of solutions =  y^(1), …, y^(N), with specific metrics such as IGD and FD relying on the known Pareto front for accuracy.

The second group of metrics assesses the quality of individual solutions y when a preference vector λ is provided.

Group 1: Metrics for Sets of Solutions

    
  * Hypervolume (HV) (↑) <cit.>: This metric evaluates both the convergence to the Pareto Front (PF) and the diversity of solutions. A low HV value indicates poor convergence, while high HV values imply better performance. The hypervolume is calculated as the volume dominated by at least one solution in the set  with respect to a reference point : 
    HV_ = Vol(  | ∃' ∈, ' ≼≼ ).

    
  * Inverted Generational Distance (IGD) <cit.>: The IGD measures the average distance between points in a reference set  and the nearest solutions in the set : 
    IGD() = 1/||( ∑_i=1^||min_y' ∈ρ(y^(i), y')^2 )^1/2.

    
  * Fill Distance (FD) <cit.>: This metric calculates the covering radius of a set of solutions , defined as the maximum minimum distance from any point in the reference set  to the nearest solution in : 
    FD() = max_y' ∈min_y ∈ρ(y, y').

    
  * Minimal Distance (l_min): This metric captures the smallest pairwise distance among all objectives:
    
    l_min = min_1 < i < j < kρ(y^(i), y^(j))

    where ρ denotes the Euclidean distance.
    
  * Smooth Minimal Distance (sl_min): This metric is a “smooth-min" version of the minimal distance function, defined as: 
    slmin = -1/h · k(k-1)log( ∑_1 < i < j < Nexp(-h ρ(y^(i), y^(j))) ).


    
  * Spacing: This metric measures the standard deviation of the minimal distances from one solution to others, with lower values indicating a more uniform distribution of objective vectors: 
    spacing = 1/N∑_i=1^N (d_i - d̅)^2,     d̅ = 1/N∑_i=1^N d_i,      d_i = min_i ≠ jρ(y^(i), y^(j)).


    
  * Span: This metric evaluates the range (span) of solutions in their minimal dimension, defined as: 
    Span = min_1 ≤ i ≤ mmax_1 ≤ k < l ≤ K | y_i^(k) - y_i^(l)|.






Group 2: Metrics for Individual Solutions

    
  * Penalty-based Intersection (PBI): This metric is a weighted sum of two distance functions d_1 and d_2, given by PBI = d_1 + μ d_2, where 
    d_1 = ⟨ y - z, λ⟩/λ,      d_2 = y - (d_1λ + z).

    
  * Inner Product (IP): This metric measures the alignment of the objective vector y with the preference vector λ: 
    IP = ⟨ y, λ⟩.

    
  * Cross Angle (ϑ): For bi-objective problems, this metric measures the angular difference between the objective vector and the preference vector: 
    ϑ = arctan(y_2 / y_1) - arctan(λ_2 / λ_1).



Those metrics are summarized in <Ref>:





 §.§ Example codes solving multitask learning problems

The following code example demonstrates solving a multiobjective fairness classification problem in four steps: defining the hyper-parameter dictionary, setting preference vectors, initializing the multitask optimizer, and solving the problem.

[style=pythonstyle, caption=An example of searching for a set of Pareto solutions.]
import torch
from libmoon.problem.mtl.core.pref_set_mtl import MTL_Pref_Solver, MTL_Set_Solver 
from libmoon.util_mtl.util import get_mtl_prefs

kwargs = 
        'batch_size' : 128,
        'lr' : 1e-3,
        'epoch' : 10,
        'solver' : 'pmgda',
        'dataset_name' : adult,
        'obj_normalization' : False,
        'n_prob' : 10,
        'sigma' : 0.6,
        'h_tol' : 1e-4,
        'device' : torch.device('cuda' if torch.cuda.is_available() else 'cpu'),
    

mtl_solver = MTL_Pref_Solver(**kwargs)
pref_mat = get_mtl_prefs(args.dataset_name, kwargs['n_prob'],
                         obj_normalization=kwargs['obj_normalization'])
res, res_history, pref_mat = mtl_solver.solve(pref_mat)




 §.§ License and usage of 
 
The license used for Adult/Compas/Credit follows Creative Commons Attribution 4.0 International (CC BY 4.0), Database Contents License (DbCL) v1.0, and CC0: Public Domain, respectively.

For academic use of , please cite our paper or GitHub. Commercial use requires author permission.   
    

unsrt  


