Profiling Pareto Front With Multi-Objective Stein Variational Gradient Descent
Year: 2022

Abstract:
Finding diverse and representative Pareto solutions from the Pareto front is a key challenge in multi-objective optimization (MOO). In this work, we propose a novel gradient-based algorithm for profiling Pareto front by using Stein variational gradient descent (SVGD). We also provide a counterpart of our method based on Langevin dynamics. Our methods iteratively update a set of points in a parallel fashion to push them towards the Pareto front using multiple gradient descent, while encouraging the diversity between the particles by using the repulsive force mechanism in SVGD, or diffusion noise in Langevin dynamics. Compared with existing gradient-based methods that require predefined preference functions, our method can work efficiently in high dimensional problems, and can obtain more diverse solutions evenly distributed in the Pareto front. Moreover, our methods are theoretically guaranteed to converge to the Pareto front. We demonstrate the effectiveness of our method, especially the SVGD algorithm, through extensive experiments, showing its superiority over existing gradient-based algorithms.

heading: 1 Introduction
text: Many scientific and engineering problems involve optimizing multiple conflicting objectives [5, 28, 4, 24], including, for example, designing wireless sensors [12], building electric power systems [30], and training neural networks with multiple tasks [35]. With multiple conflicting objectives, it is impossible to find a single solution that optimizes all the objectives simultaneously. Instead, it is essential to find a set of diverse solutions in the Pareto front that represent different preferences on the different objective functions, so that the users can have a global view of how the optimal trade-off of the different objectives look like and select the solution according to their own preference.
Unfortunately, profiling the Pareto fronts casts a key computational challenge, especially for high dimensional problems [5, 12, 28]. Traditionally, a large literature has been devoted to developing black-box, derivative-free algorithms that are suitable for black-box optimization, such as these based on evolutionary algorithms [8] and Bayesian optimization [25, 1, 39]. However, the black-box algorithms tend to be expensive and can only be applied to small scale problems due to the lack of gradient information. Gradient-based MOO algorithms have been catching attention only recently, which include mainly multiple gradient descent (MGD) based methods [10, 22, 27]. However,
∗Code is available at https://github.com/gnobitab/MultiObjectiveSampling.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
the results of these methods are still dissatisfying in finding diverse and evenly distributed Pareto solutions in complex problems.
We introduce Stein variational gradient descent (SVGD) [20, 19, 21] and Langevin dynamics [37] as efficient approaches for profiling Pareto fronts. These methods iteratively evolve a group of particles to represent the target distribution. The main difference of the two algorithms is how they distribute points. The Langevin dynamics uses stochastic noises to perturb the particle trajectories, so they can visit different areas of the Pareto front. On contrast, SVGD is a deterministic sampling algorithm that pushes the particles to high probability regions using gradient information, while enforcing diversity between the particles using a repulsive force.
In this work, we propose to a simple approach to integrate SVGD and Langevin Dynamics with MGD to draw samples from the Pareto front. Theoretical analyses are provided for both algorithms to understand their limiting distributions and their convergence speed. One challenge with MGD based sampling is that the limiting distributions do not admit explicit formulations. This is mainly because the forcing from MGD is in general not the gradient of any function. However, we can show its non-gradient component is orthogonal to a large class of functions. Assuming each objective function is strongly convex and regular, we can also show the limiting distributions concentrate on the Pareto front. Moreover, we can show the two algorithms converge to good solutions of MOO with O(1/t) and linear rate.
We test our methods on a variety of tasks, ranging from low-dimensional optimization to multi-task neural network optimization. On all the tasks tested, our method can obtain diverse and high quality Pareto solutions that distribute evenly on the Pareto front, without predefined preference vectors. Quantitatively, we substantially outperform PF-SMG and EPO with respect to the hypervolume metric.

heading: 2 Background
text: In Multi-objective optimization (MOO), we are interested in minimizing a vector-valued loss function F (x) = (f1(x), f2(x), . . . , fm(x)) ∈ Rm, where x ∈ Rd and fi : Rd → R, i ∈ [m] is the i-th scalar-valued objective function. Here we use notation [m] = {1, 2, . . . ,m}. Obviously, we can not fully optimize all the objective functions simultaneously because they may be conflicting with each other. Instead, we are interested in finding the points which can not be improved simultaneously for the objective functions, yielding the notion of Pareto optimality.
Definition 1 ( Pareto Optimality). For x1, x2 ∈ Rd, We say that x1 is dominated by x2 iff fi(x2) ≤ fi(x1),∀i ∈ [m], and F (x1) 6= F (x2). A point x∗ is called globally Pareto optimal on Rd iff it is not dominated by any other x′ ∈ Rd. A point x∗ is called locally Pareto optimal iff there exists an open neighborhood N (x∗) of x∗, such that x∗ is not dominated by any x ∈ N (x∗). The collection of globally (resp. locally) Pareto optimal points are called the global (resp. local) Pareto set. The collection of function values F (x∗) of all the Pareto points x∗ is called the Pareto front.
Our goal is to find diverse and representative solutions {xi}ni=1 from the Pareto set, so that their function values {F (xi)}ni=1 covers different preferences on different objectives. This would allow the end-users to have a global view on the optimal trade-off between the different objective functions and decide which solution based on their own preference.
Linear Scalarization One standard approach to solve MOO is using a preference vector λ = [λ1, . . . , λm] from the probability simplex on [m], i.e., S = {λ : ∑m i=1 λi = 1, λi ≥ 0, i ∈ [m]}.
Each λ ∈ S leads to a weighted objective function fλ(x) = ∑m i=1 λifi(x) and its minimizer x∗λ = arg minx fλ(x). Then as we take λ in a grid of S, we hope that the corresponding x∗λ gives a rough sketching of the Pareto front.
Although this strategy is simple and easy to implemented, it suffers from the a number of weaknesses. A key problem is that F (x∗λ) can only lie on the convex envelope of the Pareto front, and hence it only works in cases when the Pareto front is convex. In addition, a uniform grid of λ on S does not necessarily yield uniformly distributed points on the Pareto front.
Multiple Gradient Descent (MGD) MGD is a natural extension of the single-objective gradient descent to finding a Pareto point [10], which, unlike linear scalarization, works for non-convex Pareto
fronts. The idea is to iteratively update the variable x along a direction that ensures that all the objectives are decreased simultaneously (which is called Pareto improvement).
Let gi(x) = ∇fi(x) be the gradient of the i-th objective. Suppose we update the variable by x′ ← x− g∗(x), where g∗(x) is a vector field to be determined and is a small step size. By Taylor approximation, we have 〈gi, g∗〉 ≈ −(fi(x′)− fi(x))/ , which represents the decreasing rate of fi when we update x along direction g∗(x). In MGD, g∗ is chosen to maximize the slowest decreasing rate among all the objectives, that is,
g∗(x) ∝ arg max g∈Rd { min i∈[m] 〈g, gi(x)〉, s.t. ‖g‖ ≤ 1 } . (1)
Therefore, g∗(x) is encouraged to have positive inner products with all gi(x). If this is impossible to achieve, {gi(x)}mi=1 will contain conflicting directions, and we would have g∗(x) = 0 which terminates the algorithm. Using Lagrangian duality, we can show that the optimal solution to (1) is g∗(x) ∝ ∑m i=1 λ ∗ i (x)gi(x), where {λ∗i (x)}mi=1 is the solution of
min {λi} ∥∥∥∥∥ m∑ i=1 λigi(x) ∥∥∥∥∥ s.t. m∑ i=1 λi = 1, λi ≥ 0, ∀i. (2)
This optimization has a simple closed-form solution when m = 2, and a fast algorithm is offered by [35] for m > 2. By construction, when the step size is small, MGD monotonically decreases all the objectives simultaneously and will terminate when it arrives a local Pareto point; in this case we have g∗(x) = ∑m i=1 λ ∗ i (x)gi(x) = 0, suggesting that the zero vector 0 is inside the convex hull of {gi(x)}mi=1. The MGD direction g∗ provides a natural notion of gradient for multiple objectives. A point x is said to be Pareto stationary if g∗(x) = 0. Similar to the case of differentiable single-objective optimization, every Pareto local optimal is a Pareto stationary point. See [10, 35, 22] for more discussion.
However, the vanilla MGD suffers from several key weaknesses. Although it promises to return a point on the (local) Pareto set, and it is difficult to explicitly control which Pareto point it will converge to. The convergence point of MGD is implicitly determined by the initialization and the other hyper-parameters of the algorithm (e.g., step size) in a complicated way. It can be harder for MGD to converge to some Pareto points than the others. In fact, assume the Pareto front is an open set in Rd, then the vanilla MGD, when initialized from outside of Pareto front, will terminate when it reaches the boundary of the Pareto front, and hence never reach the interior points. Below is a simple example to demonstrate this.
Example 1. Consider f1(x) = x2 and f2(x) = (x − 1)2 for x ∈ R as shown in Figure 1. Then the Pareto set is the interval [0, 1]. However, when initialized outside of [0, 1], MGD will converge to either x = 0 or x = 1 and can not reach the interior points of the Pareto set unless we add a diversity-promoting mechanism, such as random noise or deterministic repulsive force. See figure 1.

heading: 3 Profiling Pareto front with Particle Dynamics
text: We improve MGD by integrating it with sampling algorithms for approximating distributions to obtain diverse solutions that cover the Pareto front more evenly. In this section, we first introduce our
main algorithm, which integrates MGD with Stein variational gradient descent (SVGD), yielding an interacting particle system that evolves diverse Pareto solutions with deterministic repulsive forces (Section 3.1). We also introduce a related idea of combining MGD with Langevin dynamics to obtain diverse solutions by injecting diffusion noises, which is found less favorable in practice than SVGD-based method but is of independent interest (Section 3.2).

heading: 3.1 Multi-Objective Stein Variational Gradient Descent
text: Stein Variational Gradient Descent (SVGD) [20] is a deterministic sampling method for finding diverse and representative points (a.k.a. particles) for approximating an un-normalized distribution, by iteratively evolving the particles with both gradient information and a special repulsive force. In this section, we integrate SVGD with MGD and adapt it to profile the Pareto front.
We derive our algorithm following how SVGD was derived in [20]. We initialize a set of points {xi}ni=1 (a.k.a. particles) and iteratively updates them to 1) mimic the MGD trajectory and 2) maintain diversity between the points. For the derivation, assume the particles {xi} follow a distribution ρ and we update them with xi ← xi − φ∗(x), ∀i ∈ [n], where φ is an update direction, selected from a predefined function spaceH, to strike the balance of the following two factors: 1) We want to choose φ to make the dynamics mimic the MGD trajectory as close as possible; this can be framed as maximizing the average inner product Ex∼ρ[〈g∗(x), φ(x)〉], to make the direction of φ and g∗ as close as possible.
2) We also want to choose φ to encourage the diversity between the particles. Assume ρ is a continuous distribution; the diversity of ρ can be measured by the entropy H(ρ), and hence we want to choose φ such that the entropy of the distribution of the updated particles, denoted by ρ′, is maximized. [20] showed that H(ρ′) = H(ρ)− Ex∼ρ[∇ · φ(x)] +O( 2), where∇· denotes the divergence operator in vector calculus. Therefore, we should choose φ to maximize−Ex∼ρ[∇·φ(x)] to encourage the diversity of the updated particles.
Overall, we optimize φ to maximize a linear combination of the two terms above: φ∗ρ = arg maxφ∈F { Ex∼ρ[〈g∗(x), φ(x)〉 −α∇ · φ(x)] } , where α is a positive coefficient that controls the importance of the divergence term. Following [20], we take F to be the unit ball of a reproducing kernel Hilbert space (RKHS) with a positive definite kernel k(x, x′). With the same derivation as [20], we get a closed form solution: φ∗ρ(x) ∝ Ex′∼ρ [g∗(x′)k(x′, x)− α∇x′k(x′, x)] .
By approximating ρ with the empirical distribution of the particles ρ = ∑ i δxi/n and iteratively updating the particles with φ∗ρ above, we obtain the following MOO-SVGD:
xi ← xi − φ̂(xi), where φ̂(xi) = 1
n n∑ j=1 g∗(xj)k(xi, xj)− α∇xjk(xi, xj). (3)
We can see that MOO-SVGD simply replaces the gradient in SVGD with the multiple gradient g∗.
The intuition of the update is clear: the first term in φ̂ makes the dynamic follow the MGD direction so that the particles are pushed towards the Pareto set when {g∗(xi)} have large magnitudes; the second term pushes the different particles away from each other, so that more areas of the Pareto set is approximated. In particular, when all the particles are in the Pareto set with g∗(xi) ≈ 0, then only the pairwise repulsive force will push the particles to be away from each other and form a uniform distribution inside the Pareto set.
If g∗ is the gradient vector of some scalar function, i.e., g∗(x) = ∇f∗(x), then our method reduces to the standard SVGD for sampling from distribution ρ∗(x) ∝ exp(−f∗(x)/α). In general, g∗ is not a gradient vector field of any scalar-valued objective function, and hence our method yields a general “non-gradient” variant SVGD. Similar to vanilla SVGD, we can characterize the evolution of the density of the particles with a nonlinear different equation. In the large particle (n → ∞) and continuous time limit, let ρt be the limit density of the particles at time t ( → 0). Following [19], the evolution of ρt is governed by,
dρt(x)
dt = ∇x · (φ∗ρ(x)ρt(x)). (4)
Unlike the case of the vanilla SVGD, the behavior of (4) is less clear due to the non-gradient nature of g∗(x). In the following, we show that the stationary distribution of (4) is connected to a Helmholtz
like decomposition of g∗(x) into a gradient component and non-gradient components. To facilitate our discussion, we define the kernel embedding of any function f with density ρ as
f[ρ](x) =
∫ k(x, y)f(y)ρ(y)dy.
Theorem 3.1. Assume ρ∗ is a fixed point of (4) for which log ρ∗(x) is continously differentiable and satisfies the Stein’s identity Ex∼ρ∗ [∇ log ρ∗(x)k(x, x′) +∇k(x, x′)] = 0 for every x′ ∈ Rd. Then we have the following decomposition of the vector field g∗:
g∗(x) = −α∇ log ρ∗(x) + (x), (5)
where ∇ log ρ∗ and  can be viewed as the gradient and non-gradient components of g∗, and (x) satisfies the following property: its kernel embedding is orthogonal to the kernel embedding of any gradient field in RKHS H, that is, let h be any C1 function with a compact support then 〈[ρ∗], ∇h[ρ∗]〉H = 0.
The existence of ρ∗ is discussed in the Appendix after the proof of Theorem 3.1. Next, we show that ρ∗ concentrates on the Pareto set. In order to do so, recall that ‖g∗(x)‖ = 0 indicates that x is on the Pareto set. Therefore we will show the norm of the MGD field g∗ is small over ρ∗ and SVGD updates can find a solution with small ‖g∗‖. To simplify our discussion, we limit our discussion on the case where a truncated (σ,md)-Gaussian kernel is used: k(x, y) = (det(2πσ2I))−
1 2 exp(− 12σ2 ‖x− y‖ 2)1‖x−y‖≤mdσ. Note that its gradient may be singular at the truncation threshold. Formally we can ignore such singularity and simply define formally: ∇yk(x, y) := x−yσ2 k(x, y)1‖x−y‖≤mdσ. We also assume the objectives functions are regular and strongly convex:
Assumption 1. Each objective function fi is C3, i.e., third-order continuously differentiable, and c1I ∇2fi c2I, |∂3j,k,lfi| ≤ c3, for i, k, l ∈ [d], where 0 < c1 ≤ c2 ≤ ∞, 0 ≤ c3 < ∞ are constants. Here ∇2fi is the Hessian matrix of fi, and A B indicates that B − A is a positive semi-definite matrix.
It is well known that, with strongly convex single objective function, gradient based algorithm can converge to the global minimum. We show that this also true for MOO-SVGD method.
Theorem 3.2. Suppose ‖g∗‖2 is C1, g∗ is L-Lipschitz and Assumption 1 holds. If we apply MOOSVGD using the (σ,md)-truncated Gaussian kernel k, the equilibrium distribution ρ∗ of MOO-SVGD satisfies
‖g∗[ρ∗]‖ 2 H =
∫ ρ∗(x)k(x, y)ρ∗(y)g∗(x)>g∗(y)dxdy
≤ ∫ ρ∗(x)k(x, y)ρ∗(y)‖g∗(x)‖2dxdy ≤ 1
2c21 (L+ α/σ2)2m2dσ 2‖1[ρ∗]‖2H.
If the MOO-SVGD density ρt following (4) is considered, for any t there is an s ≤ t such that
‖g∗[ρs](x)‖ 2 H ≤
1
2c21 (L+ α/σ2)2m2dσ 2‖1[ρs]‖ 2 H +
1 t Ex∼ρ0‖g∗(x)‖2.
Claim 1 of Theorem 3.2 says that, at the equilibrium ρ∗, the kernel embedding of g∗ is O((1 + α/σ2)2σ2) which can be arbitrarily small if σ and α are chosen to be small numbers. Claim 2 of Theorem 3.2 indicates the SVGD density flow ρt can find a solution of similar property in finite time.

heading: 3.2 Multi-Objective Langevin Dynamics
text: Besides SVGD, MGD can also be integrated with Langevin dynamics (LD), which yields diverse solutions by injecting random noise into the updates. Similar to SVGD, we can achieve this by simply replacing the gradient field in LD with g∗, yielding MOO-LD:
x← x− g∗(x) + √ 2α ξ, (6)
where ξ is a standard Gaussian noise and α is a positive coefficient that determines the magnitude of the noise. The intuition is again simple: if x is far away from the Pareto set with a large multiple gradient ‖g∗(x)‖, the dynamics will mainly drive x towards the Pareto set; when x is close to Pareto set with a small ‖g∗(x)‖, the noise term will dominate and allow x randomly diffuse on Pareto set. If g∗ is the gradient vector of some sclar-valued function f∗, that is, g∗(x) = ∇f∗(x), then its stationary distribution equals ρ∗(x) ∝ exp(−f∗(x)/α). If g∗ is not a gradient field, the Langevin dynamics is known to be irreversible [e.g., 32, 11].
As shown in our empirical results, the SVGD based method is more practically attractive than the Langevin method. We find that SVGD tends to provide more diverse and uniformly distributed points than Langevin dynamics. And for some of the challenge MOO benchmarks, Langevin dynamics can be stuck in a small local region and is slow to mix. In comparison, thanks to the deterministic updates, SVGD is found to generate mores uniformly distributed Pareto solutions in a faster speed. Another practical advantage of SVGD is that we can define the kernel to be of form k(x, y) = k0(F (x), F (y)), so that the diversity force is w.r.t. the function values {F (xi)}, rather than the variables {xi}. In the following, we consider the asymptotic property of the dynamics in (6). Consider the continuous time limit of the dynamics in (6) ( → 0), which yields a stochastic differential equation
dxt = −g∗(xt)dt+ √ 2αdwt, (7)
where wt denotes standard Wiener process and t denotes time. Assume we initialize x0 from a distribution with continuous differential density function ρ0, then the evolution the density of x(t), denoted by ρt, is governed by a Fokker-Planck equation:
d dt ρt(x) = ∇ · (g∗(x)ρt(x) + α∇ρt(x)), (8)
where∇· denotes the divergence operator in vector calculus. The stationary distribution of (7) and (8), denoted as ρ∗, uniquely exists as long as g∗(x) is coercive (e.g., if Assumption 1 holds). See e.g., [29]. It yields a decomposition of g∗(x) analogous to that of (5). Moreover, we can show the KL divergence between ρt and ρ∗ decays:
Theorem 3.3. Assume ρ∗ is a fixed point of (8) for which log ρ∗(x) is in C1(Rd). Then we have the following Helmholtz decomposition of the vector field g∗: g∗(x) = −α∇ log ρ∗(x) + (x), where  satisfies that, for any C1(Rd) function h with a compact support,  is orthogonal with∇h under ρ∗:∫ (x)>∇h(x)ρ∗(x)dx = 0. Further, if ρt(x)(x)→ 0 when ‖x‖ → ∞, the KL divergence between ρt and ρ∗ decreases monotonically with
d dt KL(ρt || ρ∗) = −αF(ρt || ρ∗), (9)
where F(ρt || ρ∗) = Ex∼ρt [‖∇ log(ρt(x)/ρ∗(x))‖ 2 ] is the Fisher divergence.
Theorem 3.3 is similar to Theorem 3.1. The differences are: 1) there is no kernel, and the decomposition is more well-understood, see for example Section 1 of [11]; 2) Interestingly, SVGD does not seem to have a similar monotonic decreasing property similar to (9). In particular, if ρ∗ is strongly log-concave, it is well known that (9) indicates that ρt converges to ρ∗ linearly in KL divergence.
Analogues to Theorem 3.2, we can also show the MOO-LD yields a distribution that concentrates on the Pareto (stationary) set. In fact, we can further show ‖g∗(x)‖ is sub-Gaussian under ρ∗.
Theorem 3.4. Under Assumption 1, let M = c23α/c1 + c22. For any constant 0 < b ≤ c14αc22 , the fixed point ρ∗ to (8) satisfies
Ex∼ρ∗‖g∗(x)‖2 ≤ 4αMd/c1, Ex∼ρ∗ exp(b‖g∗(x)‖2) ≤ 3 exp(16αbMd/c1).
The solution ρt to (8) satisfies the following
Ex∼ρt‖g∗(x)‖2 ≤ exp(− 1
2 c1t)Ex∼ρ0‖g∗(x0)‖2 + 4αMd/c1,
Ex∼ρt exp(b‖g∗(x)‖2) ≤ exp(−αbMdt)Ex∼ρ0 exp(b‖g∗(x)‖2) + 3 exp(16αbMd/c1).

heading: 4 Related Works
text: We discuss related works on MOO and also introduce the hypervolume (HV) indicator which we will use as the evaluation metric in our experiments.
Gradient-free MOO Algorithms A large volume of literature focuses on black box and derivativefree MOO, where fi(x), i ∈ {1, 2, . . . ,m} are observed through the outputs [e.g., 8, 25, 15, 39, 25, 13]. These methods are typically based on evolutionary algorithm or Bayesian optimization. However, due to the lack of gradient information, they can perform well only in low dimensional problems and can not be applied to large-scale optimization problems in modern machine learning, such as learning neural networks.
Gradient-based MOO Algorithms Another important line of work is gradient-based methods. The foundation of these methods is MGD. Representative methods include Pareto Front Stochastic Multiple Gradient (PF-SMG) [10, 22], Pareto Multi-task Learning (PMTL) [18], and Exact Pareto Optimization (EPO) [27].
PF-SMG leverages the randomness in stochastic multiple gradient descent to profile the whole Pareto front. The method keeps a set of solutions. In every optimization step, the algorithm removes dominated solutions from the set, randomly sample new solutions from the neighborhood of the survived solutions, and run MGD for another multiple steps. PF-SMG is similar to the Langevin method we discuss in Section 3.2 in that both of them leverage random noise, but PF-SMG leverages the noise in stochastic gradient, which will vanish to zero when the step size decreases to zero. Therefore, PF-SMG would reduce to MGD except the selection step. In constrast, the diffusion noise converges to Brownian motion with small step size.
PMTL and EPO are designed to find local Pareto optimal solutions that satisfy specific user preference. They can also be used to trace the Pareto front by varying the user preference vector uniformly. However, sampling the preference vector uniformly will not necessarily result in a uniform solution set on the Pareto front, which is what we observe in the experiments. Compared with PTML and EPO, our method does not require pre-defined user preference and can obtain more diverse and uniformly distributed points as we show in experiments.
Pareto Hypernet Methods Several methods, such as [31, 17, 33], proposed to learn a neural network, called Pareto hypernet, which can generate pareto-optimal solutions given a preference vector as input. In comparison, our methods directly generates a set of solutions (or particles) that covers the Pareto front. The usage of neural network can be helpful in some cases, but it also makes the algorithm more complicated and less transparent, and the choice of neural network structure has heavy influence on the result. In comparison, our algorithm is much simpler, transparent and directly outputs the solutions that we want. More empirical evidence and discussion can be found in the Appendix.
Hypervolume The hypervolume (HV) indicator is a standard metric for evaluating the quality of sets of solutions in MOO. Assume r = [r1, . . . , rm] is an reference point that is an upper bound of the objectives, such that supx fi(x) ≤ ri, ∀i ∈ [m]. For a given set of solutions X = {xi}ni=1, its hypervolume indicator HV (X ) is the measure of the region between F (X ) and r,
HV (X ) = Λ ({ q ∈ Rd | ∃x ∈ X : q ∈
m∏ i=1 [fi(x), ri]
}) , (10)
where Λ(·) denotes the Lebesgue measure. Directly optimizing the hypervolume metric can be problematic, because it is piece-wise constant. Hypervolume indicator gradient ascent (HIGA) [36] relaxes the hypervolume indicator to get usable gradient for optimization, and a recent work [9] is mathematically similar to HIGA. However, [36] and [9] have two main drawbacks compared with our sampling-based method: (1) its results depend on the choice of reference point; (2) it may get stuck in bad local minima. See the Appendix for more empirical comparison.

heading: 5 Experiments
text: We show our method, Multi-Objective Stein Variational Gradient Descent (MOO-SVGD) can efficiently profile the Pareto front on a variety of problems, including the ZDT problem set, tri-objective problem MaF1, trade-off between accuracy and fairness metrics in fair ML, as well as multi-task
learning on MultiMNIST, MultiFashion and MultiFashion+MNIST. Our method outperforms EPO and PF-SMG, two previous baselines for profiling the Pareto front with gradients. We also include Multi-Objective Langevin Dynamics (MOO-LD) for comparison, and find that MOO-SVGD yields much better empirical performance against MOO-LD. MOO-LD in general needs longer iterations to cover the Pareto Front. Please refer to the Appendix for more discussion. We leave the detailed algorithm configuration in the Appendix due to space limitation.
ZDT Bi-objective Problems We compare our methods with other baseline methods on the ZDT problem set [40]. The ZDT problems have a 30-dimensional variable and two objective functions. The details of ZDT problems can be found in the Appendix. In ZDT problems, both objective functions are differentiable. We choose ZDT-1, ZDT-2, and ZDT-3 to demonstrate the performance of our algorithm when the Pareto front is convex, concave, and discontinuous. For baselines, we compare with naive optimization on hypervolume , state-of-the-art gradient-based MOO algorithms (EPO [27] PF-SMG [22]) and Pareto hypernet algorithms (COSMOS [33]). All the results of the baselines are reproduced using their official implementation. For fair comparison between black-box and white-box optimization algorithms, we report the increase of the hypervolume indicator (Eq. (10)) over the wall-clock time. For fairness of comparison, we use the same reference point when computing hypervolume indicator for every algorithm. Moreover, to make the comparison of wall-clock time fair, we run all the algorithms with CPU on the same computational platform with 192GB memory and 48 cores. The results are shown in Fig. 2. We constrain the maximal wall-clock time in the figure to be 60 seconds. We use 50 particles for MOO-SVGD. Because EPO and COSMOS need user preference vectors, we uniformly sample 50 unit vectors from a semi-circle to trace a Pareto front, following the strategy provided by the authors [27, 33]. We use a two-layer ReLU network with 1,000 hidden units for the hypernet in COSMOS.
As the figures shows, our method, MOO-SVGD, obtains the largest hyper-volume in all three problems. In the limit of 60 seconds, EPO generally has the best performance among all the three baselines. Naive optimization of hyper-volume get similar performance as the stochastic method, PF-SMG, though the latter is much slower than the former. COSMOS cannot cover the entire Pareto
front. In ZDT1 and ZDT3, MOO-SVGD converged to a better solution. In comparison, another sampling based method, MOO-LD, can only cover a limited subset of the whole Pareto front.
MaF1 Tri-objective Problem We test our algorithm on MaF1 [3], a problem with three objectives to showcase its ability when there are more than two objectives. The objective functions are defined on [0, 1]d. The specific expression of f1, f2 and f3 can be found in Appendix due to limited space. The Pareto front of the problem can be analytically solved, which is {x : f1(x)+f2(x)+f3(x) = 2}. The reference point for calculating hypervolume is (1.0, 1.0, 2.0). The results are shown in Fig 3 and Table. 1. All the three methods, EPO, MOO-LD, and MOO-SVGD can find Pareto optimal solutions. However, MOO-LD can only cover a very small subset of the Pareto front. EPO can find diverse Pareto optimal solutions, but they are not distributed uniformly. MOO-SVGD obtains the most visually uniform solution set, and hence reaches the highest performance w.r.t. hypervolume. MOO-LD is not tested in the following parts due to its poor performance on toy experiments.
Trade-off Between Accuracy and Fairness When a ML system is applied to a real-world scenario, it may be biased against people with sensitive attributes, e.g., their genders or races. This leads to unfairness in ML, which is a key issue in real applications of ML. A common practice to avoid unfairness is to add fairness regularization during training. Unfortunately, introducing fairness regularization may hurt the prediction accuracy of the model. The user may have various requirements for the trade-off between accuracy and fairness in different environments. This can be framed into a MOO problem, which simultaneously optimize the typical training loss and the fairness penalties.
For this experiment, please refer to the appendix to see how we construct our dataset. We choose gender to be the sensitive attribute, and we use disparate impact as the fairness criterion. Disparate impact can be measured by the CV score [2], CV(f) = |p(ŷ = 1|A = 0)− p(ŷ = 0|A = 0)|, where f is our classifier, ŷ is the prediction, and A is the sensitive attribute. ŷ = 1 if f(x) > 0; otherwise, ŷ = 0. CV score implies the gap between the probabilities of getting positive outcomes in different sensitive groups. As in [23], we use a linear classifier with 87 parameters. The classification loss `1 is the typical binary cross-entropy loss. For fairness penalty, instead of directly optimizing the CV score, we use a convex surrogate of CV score provided in [23], `2(f) = [Exi,ai∈D(ai − ā)f(xi)]
2, where xi is the feature of instance i; ai is the sensitive attribute of instance i; D is the dataset, and ā is the averaged value of all ai in the dataset. Basically, it is minimizing the empirical covariance between the the prediction and the sensitive attribute.
We compare with PF-SMG [23] and EPO [27]. The number of particles is 10 for MOO-SVGD and EPO. For PF-SMG, we use two different versions. The first version has a maximal particle number of 600; the second is 50. The two versions are denoted by PF-AMG-50 and PF-SMG-600. The results are shown in Fig. 4. Generally, MOO-SVGD yields the best solution set at convergence. Among all
the methods, only MOO-SVGD and PF-SMG-600 can get a solution set that is uniformly distributed on different preferences. PF-SMG-50, due to poor exploration ability of the random strategy, shrinks to the two end points of PF-SMG-600. EPO totally fails in this task. Uniformly sampling preference vector from the semi-circle is problematic, since the Pareto front lies in a very small region of the semi-circle, and most of the preference vectors do not ends up with an ‘exact’ solution at all.
Multi-task Learning with Neural Networks MultiMNIST, MultiFashion and MultiFashion+MNIST, introduced in [34], are three standard benchmarks for multi-task deep learning. The datasets are created in a similar way: each image in the datasets is a combination of two different images. One of them is on the top-left, and the other is on the bottom-right. In MultiMNIST, both images are from the MNIST dataset [16]. Similarly, MultiFashion uses two images from the FashionMNIST dataset [38]. For MultiFashion+MNIST, the top-left image is from MNIST, while the bottom-right one is from FashionMNIST. All three datasets have 120,000 samples in the training set, and 20,000 samples for test set. These datasets are used in [18, 27, 26].
The baseline methods are linear scalarization, EPO [27] and Pareto Multi Task Learning (PMTL) [18]. EPO and PMTL both require preference vectors to work. We follow the setting in [27], where the authors use LeNet [16] and 5 preference vectors uniformly sampled from a semi-circle. We also use 5 particles for MOO-SVGD. We find that training with MOO-SVGD can be problematic in the early stage, where networks with different initialization can be too close to each other and cause gradient explosion. Therefore, we initialize the networks with same parameters, then warm-up them by training with naive linear scalarization objective for one epoch, where the weights are sampled uniformly from [0, 1]. We switch back to MOO-SVGD in the remaining epochs. For all the experiments, all the methods are trained for 100 epochs with a batch size of 256. Moreover, they are optimized with SGD and the learning rate is 0.001.

heading: 6 Conclusion, Limits, Impacts, and Open Questions
text: In this paper, we propose Multi-objective Stein Variational Gradient Descent (MOO-SVGD). Empirically, it can trace the Pareto front evenly without user preference. There are a number of limits that deserve future efforts. On the practical side, since our methods leverage SVGD and Langevin dynamics, they necessarily inherent their limits, such as the sensitivity on the choice of kernel and step sizes; studies on automatizing the selection of these hyper-parameters can improve the method significantly. Theoretically, analyzing our methods are challenging due to the lack of explicit form of the invariant distributions. For example, although we show that the limiting distributions concentrate on the Pareto front, it is an open question to theoretically characterize in what sense they are evenly distributed in the Pareto front. Moreover, it would be interesting to draw theoretical understandings on why the MOO-SVGD performs better than the MOO-Langevin dynamics.
In terms of social impact, as a new computational tool for MOO, it allows us to better trade-off different objectives in machine learning and hence improve ML systems in both accuracy and various trustworthy metrics such as fairness, safety, as we demonstrate with our experiments. We will release code online to promote the applications of our approach, although it does open the possibility of malicious use of our techniques for adversarial purposes.

heading: Acknowledgements
text: The work is conducted in the statistical learning and AI group in computer science at UT Austin, which is supported in part by CAREER-1846421, SenSE-2037267, EAGER-2041327, and Office of Navy Research, and NSF AI Institute for Foundations of Machine Learning (IFML). Xingchao Liu is supported in part by a funding from BP. X. T. Tong’s research is supported by Singapore MOE Academic Research Funds R-146-000-292-114. The authors thank the reviewers for all the suggestions made in the reviewing process.

Review summaries:
This paper proposes a method to learn multiple Pareto solutions for multi-objective optimization (MOO). Specifically, the method iteratively updates a set of points by Stein variational gradient descent (SVGD) to encourage the diversity between the solutions. Another counterpart method based on Langevin dynamics is also provided. The methods are theoretically guaranteed to converge to the Pareto front. The experimental results demonstrate the superiority of the proposed SVGD type of method compared with EPO. The superiority is shown visually and by the hypervolume of the solutions.

In this work the authors propose a gradient-based method for identifying the Pareto Front of a Multi-Objective Optimization problem.   More specifically, they propose two methods, one based on Stein Variational Gradient Descent, and another based on overdamped Langevin dynamics to evolve forward a set of particles to the (local) pareto front.  The approach is based on the multiple gradient descent (MGD) formulation where the next particle descent direction is chosen to maximise the slowest decreasing rate amongst the objectives.  Assuming particles are advected along an unknown vector field $\phi^*$ the authors introduce a variational formulation which selects the $\phi$ according to MGD objective while encouraging diversity through the usual entropy regularisation.   This gives rise to SVGD-like dynamics.   An analogous formulation for Langevin Dynamics is also proposed.  The authors provide some theoretical guarantees for both in terms of concentration of the dynamics to the Pareto front.  Results are demonstrated for a number of numerical experiments. 

This paper studies how to acquire diverse and representative Pareto solutions by gradient-based methods. Two algorithms -- MOO-SVGD and MOO-LD are proposed by using stein variational gradient descent and Langevin dynamics respectively, with convergence (to local Pareto optimal) guarantees. The experimental results demonstrate that MOO-SVGD can obtain visually uniform solution sets of ZDT and MaF1 problems, and get better Pareto fronts on fair ML and multi-task learning.