# Stochastic Tangential Pareto Dynamics Provably Samples the Whole Pareto Set

Zachary Jones, Pietro Marco Congedo, Olivier Le Maitre

► To cite this version:

Zachary Jones, Pietro Marco Congedo, Olivier Le Maitre. Stochastic Tangential Pareto Dynamics Provably Samples the Whole Pareto Set. 2025. hal-04954830

HAL Id: hal-04954830

https://inria.hal.science/hal-04954830v1

Preprint submitted on 18 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 STOCHASTIC TANGENTIAL PARETO DYNAMICS PROVABLY

2 SAMPLES THE WHOLE PARETO SET *

3 ZACHARY JONES, PIETRO MARCO CONGEDO, OLIVIER LE MAÎTRE

4
5
6
7
8
9
10
11
12
13
14
15

**Abstract.** The framework of stochastic multi-objective programming allows for the inclusion of uncertainties in multi-objective optimization problems at the cost of transforming the set of objectives into a set of expectations of random quantities. The stochastic multigradient descent algorithm (SMGDA) gives a solution to these types of problems using only noisy gradient information. However, a bias in the algorithm causes it to converge to only a subset of the whole Pareto front, limiting its use. We analyze the source of this bias and prove the convergence of SMGDA to a stationary point in the nonconvex L-lipschitz smooth case. First, based on this analysis, we propose to reduce the bias of the stochastic multi-gradient calculation using an exponential smoothing technique. We then propose a novel approach to exploring the whole Pareto set by combining the unbiased stochastic multigradient with an additive non-vanishing noise that guides the dynamics of the iterates tangential to the Pareto set. We finish by proving that our algorithm, Stochastic Tangential Pareto Dynamics (STPD), generates samples concentrated on the whole Pareto set.

16
17
18
19
20
21
22
23
24
25
26
27
28

**Key words.** Stochastic Programming, Multi-Objective Optimization, Stochastic Multi-Objective Optimization, Gradient, Sampling, Pareto Front, Pareto Set

18 MSC codes. 68Q25, 68R10, 68U05

1. **Introduction.** Finding all possible solutions to an optimization problem with competing objectives and uncertainty requires specialized methods. In general there is no unique solution to an optimization problem with competing objectives but rather a collection of possible tradeoffs. A Pareto optimal point is then a design point which cannot be improved along any objective without degrading the performance of others. *the Pareto set* is the set of all Pareto optimal points. Introducing stochasticity into the objective functions complicates the process of finding the Pareto set. The objective functions then become random variables making both optimization and the determination of Pareto optimality into probabilistic tasks.

28
To handle the challenges of stochastic multi-objective optimization we work within

<footnote>*Submitted to the editors DATE.</footnote>

<footnote>**Funding:** This work was funded by Horizon NEXTAIR.</footnote>
the framework of stochastic programming. We formulate our stochastic multi-objective optimization problem as the minimization of a set of expected values of random quantities of interest. This changes the stochastic problem into a deterministic one. Unfortunately, we rarely have access to analytic expressions for the expected values of the quantities of interest. Furthermore, estimating them can prove to be computationally costly. Therefore, we turn to stochastic approximation.

The stochastic multi-gradient descent algorithm (SMGDA) is the extension of stochastic approximation methods to the multi-objective context. It has attracted attention in engineering, operations research, and machine learning fields [5, 11, 9, 4]. At each iteration SMGDA uses stochastic gradient information from each objective to compute a direction of descent, the stochastic multigradient. This approach has two main benefits. SMGDA iterations can be performed with as few as a single gradient sample from each objective, making estimation unnecessary. It also has the advantage of provably converging to the Pareto set [9]. However, the stochastic multigradient has a known bias [5, 11]. Because of this bias, previous analysis of the convergence of SMGDA beyond the strongly convex setting have depended on strong assumptions or even concluded that the algorithm does not converge [11, 5, 4].

Our first contribution is to analyze the source and effect of this bias. We prove that, despite bias in the stochastic multigradient, SMGDA converges to a *subset* of the Pareto set in the convex L-Lipschitz smooth setting. We additionally show convergence of SMGDA in the strongly convex and nonconvex L-Lipschitz smooth settings.

Our second contribution is a sampling approach to generate candidate Pareto optimal points, stochastic tangential Pareto dynamics (STPD). It has two main ingredients, a debiased stochastic multigradient step and an additive noise term. To bias the stochastic multigradient, based on our analysis of SMGDA, we propose a method using exponential smoothing. We then combine this debiased stochastic multigradient with a nonvanishing noise term which adds noise perpendicular to the expected stochastic multigradient. This noise term guides the sampling procedure tangentially along the Pareto set. We go on to prove that STPD generates samples
which lie arbitrarily close to the whole Pareto set in finite time.

Our companion paper then tackles the complementary problem of probabilistic inference of the Pareto set using a noisy fixed-size sample set [CITATION NEEDED].

We begin the manuscript by recalling some basic concepts in both deterministic and stochastic multi-objective optimization in Section 2. We then introduce the scalar potential for multi-objective optimization in Section 3.1, and use it to prove the convergence of MGDA in Section 3.2. We use the foundation built in Section 3 to make new, more general, proof of convergence of SMGDA in Section 4.1. We also see that SMGDA converges only to a subset of the Pareto front. We examine the bias of SMGD in Section 4.2. We introduce our sampling approach, Stochastic Tangential Pareto Dynamics, in Section 5. We continue to show that it samples the whole of the Pareto front under mild assumptions and offer discussion about the potential applications, pitfalls, improvements, and variations to our proposed approach. Finally, we give our conclusions and future outlook in Section 6

## 2. Background on Gradient Based (Stochastic) Multi-Objective Optimization.
We start by introducing stochastic programming before reviewing the multigradient descent algorithm and Pareto concepts.

Let us look at the stochastic multi-objective programming problem. Given a probability space $(\Theta, \mathcal{F}, \mu)$ and a set of $k$ quantities of interest collected as a vector valued function $F := [f_i(\mathbf{x}, \mathbf{W}(\theta)), \dots, f_k(\mathbf{x}, \mathbf{W}(\theta))] : \mathcal{X} \subseteq \mathbb{R}^d \times \mathbf{W}(\Theta) \mapsto \mathcal{Y} \subseteq \mathbb{R}^k$, our goal is to find all $\mathbf{x}^*$ such that

$$ (2.1) \qquad \mathbf{x}^* \in \underset{\mathbf{x} \in \mathcal{X}}{\operatorname{argmin}} \{\mathbb{E}[f_i(\mathbf{x}, \mathbf{W}(\theta))], \dots, \mathbb{E}[f_k(\mathbf{x}, \mathbf{W}(\theta))]\} $$

This formulation is sufficiently general to incorporate mean/variance minimization, conditional value at risk, and quantile functions as objectives. Notably, the problem is deterministic. If we have analytical expressions for $G(\mathbf{x}) = [g_1(\mathbf{x}), \dots, g_k(\mathbf{x})]$, $\mathcal{X} \mapsto \mathcal{Y} \subseteq \mathbb{R}^k$ with $g_i(\mathbf{x}) = \mathbb{E}[f_i(\mathbf{x}, \mathbf{W}(\theta))]$, we can treat this as a deterministic problem
and solve it as

$$ (2.2) \qquad \mathbf{x}^* \in \operatorname*{argmin}_{\mathbf{x} \in \mathcal{X}} \{g_1(\mathbf{x}), \dots, g_k(\mathbf{x})\}. $$

Clearly we will rarely have a unique solution to problem 2.1. Instead, given two solutions, $\mathbf{y}', \mathbf{y} \in \mathcal{Y}$, $\mathbf{y}'$ dominates $\mathbf{y}$ ($\mathbf{y}' < \mathbf{y}$) if it is lower componentwise in all dimensions, e.g. $\forall j \in \{1, \dots, k\} \mathbf{y}'_i < \mathbf{y}_i$. If no $\mathbf{y}'$ exists which dominates a point $\mathbf{y}$, we say it is *undominated*. The Pareto front is defined as the set of undominated points.

DEFINITION 2.1 (Pareto Front).

$$ (2.3) \qquad P(\mathcal{Y}) := \{\mathbf{y} \in \mathcal{Y} \mid \exists \mathbf{y}' \in \mathcal{Y} \setminus \mathbf{y} \text{ s.t. } \mathbf{y}'_j \le \mathbf{y}_j \forall j \in \{1, \dots, k\}\}. $$

This set gives a mathematical description of Pareto optimality. If $\mathbf{y} \in P$ it means that, in at least one pair of criteria, there is no global performance improvement. We can define the set of solutions in design space, the Pareto set, as the pre-image of the Pareto front of $G$.

DEFINITION 2.2 (Pareto Set).

$$ (2.4) \qquad P^{-1}(\mathcal{X}) := (P \circ G)(\mathcal{X}) = \{\mathbf{x} \in \mathcal{X} \mid G(\mathbf{x}) \in P(\mathcal{Y})\}. $$

Points $\mathbf{x}^*$ in 2.2 can be found using the multigradient descent algorithm.

It is an extension of gradient descent to the multi-objective context using the multigradient.

At each iteration we calculate

$$ (2.5) \qquad \nabla_x^C \{G\}(\mathbf{x}_t) = \underset{\mathbf{d}}{\operatorname{argmax}} \max_i \langle \nabla g_i(\mathbf{x}_t), \mathbf{d} \rangle + \frac{1}{2} \|\mathbf{d}\|_2^2. $$

Under sufficient conditions on $G$, the sequence of points $\mathbf{x}_{t+1} = \mathbf{x}_t - \varepsilon \nabla_x^C \{G\}(\mathbf{x}_t)$ will then converge to a Pareto stationary point
106 DEFINITION 2.3 (Pareto Stationary Point). $\mathbf{x}^*$ such that

$$ (2.6) \qquad \nabla_x^C \{G\}(\mathbf{x}^*) = \mathbf{0} $$

108 In the case where all objective functions are convex this forms a sufficient condition
109 for Pareto optimality [6].

110 In later sections, we will use Pareto stationarity as a tool to search for Pareto
111 optimal points.

112 In practice, solving Eq. 2.5 in its primal form is cumbersome, and we will use its
113 dual formulation, which reduces to a quadratic subproblem in the number of objec-
114 tives,

$$ (2.7) \qquad \nabla_x^C \{G\}(\mathbf{x}) = \sum_i \alpha_i^* \nabla g_i(\mathbf{x}) \\ \text{s.t.} \quad \alpha^* = \underset{\alpha \in \Delta^{k-1}}{\operatorname{argmin}} \left\| \sum_i \alpha_i \nabla g_i(\mathbf{x}) \right\|_2^2, $$

118 where $\Delta^{k-1}$ is the $k-1$ dimensional simplex. This subproblem has a closed form
119 solution in the bi-objective case and can be efficiently solved for more objectives
120 using the Franke-Wolfe algorithm [3]. The multigradient, $\nabla_x^C \{G\}(\mathbf{x})$, can be concisely
121 understood as the minimum-norm convex combination of gradients, as seen in Fig. 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Examples of Calculation of Direction of Descent $\nabla_x^C \{G\}(\mathbf{x})$ for two objectives.</figcaption></figure>

122

123 3. Solving Deterministic Multi-Objective Problems with MGDA. In
124 this section we reformulate and solve the deterministic multi-objective problem. To
reformulate the multi-objective problem we make an equivalent problem with only one
objective, a scalar potential. This provides a new way prove convergence of MGDA
and to view multi-objective optimization. In Section 4 we will extend these techniques
to the stochastic case and in Section 5 we will exploit this framework to prove that
our proposed approach samples points along the whole Pareto set.

3.1. Potential function for Multi-Objective Optimization. To reformu-
late the multi-objective problem we make an equivalent problem with only one objec-
tive, a scalar potential. Previous attempts have used the coefficients, α*, defined in
Eq. 2.7 to define a pseudo-objective Ψ(x) = Σᵢ αᵢ* gᵢ(x) [5]. However, the parameters
α* are not constant with respect to x and so the gradient, ∇ₓ_c{G}(x), is not ∇Ψ(x).
Consequently Ψ(x) is minimized only at the minima of each individual objective and
not over the whole front.

We propose a potential function that has minima over the whole of the Pareto
set. We define the potential, $\Phi_G(\mathbf{x}) : \mathbb{R}^d \to \mathbb{R}$, as follows:

139 DEFINITION 3.1 (Equivalent Scalar Potential Function).

$$
(3.1) \quad \Phi_G(\mathbf{x}) = \int_0^1 \langle \nabla_x^C \{G\}(s\mathbf{x}), \mathbf{x} \rangle ds = \int_0^1 \sum_i \alpha_i^* (s\mathbf{x}) \langle \nabla g_i(s\mathbf{x}), \mathbf{x} \rangle ds.
$$

141 The potential Φ<sub>G</sub>(x) is C<sup>1</sup> smooth, locally Lipschitz continuous, and it can be seen
142 that ∇Φ(x) = ∇<sub>x</sub><sup>C</sup>{G}(x) by the gradient theorem. Clearly this function has minima
143 at Pareto stationary points.

144 It can be readily seen in Fig. 2d that the minima of Φ<sub>G</sub> lies on the Pareto set.
145 In the case where there are nonconvex objectives, the set of Pareto stationary points
146 can form a saddle point in Φ<sub>G</sub>(x). The global minima remains on the Pareto set. We
147 remark here that we have used the curve from 0 to x parameterized by t to perform
148 the line integral that defines the potential function. This is done for convenience only.
149 Any continuous curve in X with endpoint x can be used so long as the initial point
150 is consistent and X is pathwise connected.

151 Using 3.1, we can make claims about the set of objective functions using a single
152 scalar value.
<figure><img src="image_3.png" /><figcaption>Fig. 2: Top row: Competing objective functions in the cases with local minima (left) and bi-convex objectives(right). Bottom row: Minima of the equivalent scalar potential function lie on the Pareto front, with Pareto stationary nonoptimal points revealing themselves as saddle points (left).</figcaption></figure>

**3.2. Convergence of MGDA.** We use the scalar potential to make a new proof of the convergence of MGDA. In the next section, we will use a similar strategy to show the conditions required for SMGDA to converge.

The minima of $\Phi(\mathbf{x})$ are the set of Pareto stationary points, and performing gradient descent on $\Phi(\mathbf{x})$ yields Pareto stationary points.
THEOREM 3.2. Let $\Phi_G(\mathbf{x})$ be a potential function defined as in 3.1. With $\Phi_G$:

$\mathcal{X} \mapsto \mathcal{Y}$, $\mathcal{X}$ compact. Then the sequence

$$ (3.2) \qquad \mathbf{x}_{t+1} = \mathbf{x}_t - \varepsilon \nabla \Phi_G(\mathbf{x}). $$

converges to a Pareto stationary point.

Proof. From the assumption that $\Phi_G$ is L-Liptshitz smooth

$$ (3.3) \qquad \Phi_G(\mathbf{x}_{t+1}) \le \Phi_G(\mathbf{x}_t) + \langle \nabla \Phi_G(\mathbf{x}_t), \mathbf{x}_{t+1} - \mathbf{x}_t \rangle + \frac{L}{2} \|\mathbf{x}_{t+1} - \mathbf{x}_t\|^2 $$

$$ (3.4) \qquad \le \Phi_G(\mathbf{x}_t) - \varepsilon\left(1 - \frac{L\varepsilon}{2}\right) \|\nabla_x^c\{G\}(\mathbf{x}_t)\|^2. $$

Picking $\varepsilon \le \frac{2}{L}$, summing, using the fact that $\nabla\Phi(\mathbf{x}) = \nabla_x^c\{G\}(\mathbf{x})$, and dividing by $T$ gives us the result

$$ (3.5) \qquad \frac{1}{T} \sum_{t=1}^{T} \|\nabla_x^c\{G\}(\mathbf{x}_t)\|^2 \le \frac{\Phi_G(\mathbf{x}_0) - \Phi_G(\mathbf{x}_{T+1})}{T\varepsilon}. \quad \square $$

4. Solving Stochastic Multi-Objective Problems with SMGDA. We now turn our attention to the stochastic formulation of the multigradient descent algorithm, SMGDA. It has already been proven that SMGDA converges to the Pareto front almost surely under strong assumptions. Previous approaches either require the individual objective functions are almost surely convex with at least one being almost surely strongly convex or they assume that the bias in the stochastic multigradient vanishes in the long time limit.[9, 4, 5]. These cases are not realistic. It has also been shown that there exist cases where there is a nonvanishing bias in the stochastic multigradient along the whole Pareto set [11]. The authors went on to conclude that SMGDA does not converge in these cases. However, this is not necessarily true.

Here, we extend the tools we have developed in the previous section to the stochastic case. This allows us to give proof of convergence of SMGDA in the nonconvex and biased case. We then go on to show that, when Pareto stationary points for which the stochastic gradient estimator is unbiased exist, that SMGDA converges to one of
those points.

In the stochastic programming setting the set of objective functions {$g_i$}$_{i=1,...,k}$ now have the functional form of an expectation $g_i(\mathbf{x}) = \mathbb{E}[f_i(\mathbf{x}, \mathbf{W}(\theta))]$. Unfortunately, we do not often have analytic expressions for our objectives and must work with realizations of the quantities of interest, {$f_i(\mathbf{x}, \mathbf{W}(\theta))$}$_{1,...,k}$ and their gradients. It is possible to estimate the value of our objectives and plug them in to a deterministic optimization algorithm. However, this strategy introduces residual error and can require costly resampling at each iteration. Ignoring the random nature of the problem by naively substituting samples of the quantities of interest, $f_i(\mathbf{x}_t, \mathbf{W}(\theta)_t)$, for their true value in a deterministic multi-objective method could lead to convergence to a suboptimal point.

Instead, we estimate the gradient of each objective, $\nabla g_i(\mathbf{x})$, using a single realization of the gradient of the corresponding quantity of interest $\nabla f_i(\mathbf{x}, \mathbf{W}(\theta))$. These estimates are then plugged into Eq. 2.7 to calculate the stochastic multigradient $\nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta))$. It is assumed that $\nabla f_i(\mathbf{x}, \mathbf{W}(\theta))$ is an unbiased estimator of $\nabla g(\mathbf{x})$ with finite second moment, e.g. $\mathbb{E}[\nabla f(\mathbf{x}, \mathbf{W}(\theta))] = \nabla g(\mathbf{x})$ and $\mathbb{V}(\nabla f(\mathbf{x}, \mathbf{W}(\theta))) < \infty$ with $\mathbb{V}(\cdot)$ as the variance functional. The heart of the stochastic multigradient descent algorithm is then the sequence of points generated by the relation

$$ (4.1) \qquad \mathbf{x}_{t+1} = \mathbf{x}_t - \varepsilon_t \nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}). $$

We require that $\varepsilon_t$ meet the requirement that

$$ (4.2) \qquad \frac{\sum_t^T \varepsilon_t^2}{\sum_t^T \varepsilon_t} \xrightarrow{T \to \infty} 0. $$

This is a standard requirement in stochastic optimization [1, 9].

It is straightforward then to extend the potential function from Section 3 to the stochastic multi-objective case.
207

**DEFINITION 4.1 (Stochastic Equivalent Scalar Potential Function).**

$$
208
(4.3) \qquad \Phi_F(\mathbf{x}) = \int_0^1 \langle \mathbb{E}[\nabla_x^C\{F\}(s\mathbf{x})], \mathbf{x} \rangle ds.
$$

209
210

The gradient of this potential function can be approximated using individual realiza-
tions of **W**(**θ**) as **∇**<sub>**x**</sub><sup>C</sup>{**F**}(**x**, **W**(**θ**)).

These quantities are summarized in Tab. 1.

<table><caption>Table 1: Quantities used in MGDA and their stochastic programming equivalents in SMGDA.</caption>
  <thead>
    <tr>
      <th></th>
      <th>Objectives</th>
      <th>Gradients</th>
      <th>Potential</th>
      <th>Direction of Ascent</th>
      <th>step size</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>MGDA</td>
      <td>g(<b>x</b>)</td>
      <td>∇g(<b>x</b>)</td>
      <td>Φ<sub>G</sub></td>
      <td>∇<sub><b>x</b></sub><sup>C</sup>{G}(<b>x</b>)</td>
      <td>ε</td>
    </tr>
    <tr>
      <td>SMGDA</td>
      <td>E[f(<b>x</b>, <b>W</b>(θ))]</td>
      <td>∇f(<b>x</b>, <b>W</b>(θ)<sub>i</sub>)</td>
      <td>Φ<sub>F</sub></td>
      <td>∇<sub><b>x</b></sub><sup>C</sup>{F}(<b>x</b>, <b>W</b>(θ))</td>
      <td>ε<sub>t</sub></td>
    </tr>
  </tbody>
</table>

211

212
213
214
215
216

**4.1. Convergence of SMGDA.** We can now show new proof of the conver-
gence of the SMGD algorithm to a Pareto stationary point in the nonconvex case
with biased gradients. Proofs of convergence for convex and strongly convex ob-
jective functions in the biased setting with bounded gradients can be found in the
appendix.

217

We first make an assumption on the bias of the gradient estimator,

218

ASSUMPTION 4.2. ∃a > 0, b ≥ 0 such that

$$
219
(4.4) \quad \langle \mathbb{E}[\nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta))], \nabla_x^C\{G\}(\mathbf{x}) \rangle \ge a \|\nabla_x^C\{G\}(\mathbf{x})\|_2^2 - b.
$$

220

*almost everywhere.*

221
222
223
224
225
226
227

This assumption intuitively states that the stochastic gradient is "in-line" with the
true gradient up to a constant bound, *b*. Many general proofs of the convergence of
SGD in the single objective case make a similar assumption, setting *b* = 0 directly [1].
This assumption is more general. We can always find *b* ≥ 0 such that this relationship
holds with arbitrarily high probability. As we will see later, if we can find a pair *a*, *b*
such that *b* = 0 and Assumption 4.2 holds almost surely over all *x* ∈ 𝒳, then SMGD
converges to a Pareto stationary point. The interested reader can pursue more details
228 in the appendix.

229 We also assume that

230 ASSUMPTION 4.3. $\exists a_G > 2a - 1$ and $b_G \ge 0$ such that

$$231 \quad (4.5) \quad \mathbb{E}[||\nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta))||_2^2] \le a_G ||\nabla_x^C\{G\}(\mathbf{x})||_2^2 + b_G$$

232 almost everywhere.

233 Assumptions 4.2 and 4.3 are not strong and imply two things. The first, that the effect of the bias is lower bounded:

$$235 \quad (4.6) \quad \mathbb{E}[\langle \nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta)), \nabla_x^C\{G\}(\mathbf{x}) \rangle - ||\nabla_x^C\{G\}(\mathbf{x})||_2^2] \ge (a-1) ||\nabla_x^C\{G\}(\mathbf{x})||_2^2 - b.$$

236 Second, setting $M_V \ge (a_g + 1 - 2a)$ and $M_0 \ge 2b + b_G$, we see that assumptions 4.2 and 4.3 imply 4.7.

$$238 \quad (4.7) \quad \mathbb{E}[||\nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta)) - \nabla_x^C\{G\}(\mathbf{x})||_2^2] \le M_V ||\nabla_x^C\{G\}(\mathbf{x})||_2^2 + M_0$$

239 Readers may recognize Eq. 4.7 as a common assumption in stochastic approximation [1]. We now prove the convergence of SMGDA to a Pareto stationary point in the nonconvex setting.

242 THEOREM 4.4. Let $F = \{f_i(\mathbf{x}, \mathbf{W}(\theta))\}_{i=1,...,k}$ be a vector valued function $F$ : $\mathcal{X} \to \mathcal{Y} \subseteq \mathbb{R}^k$ with each $g_i(\mathbf{x}) = \mathbb{E}[f_i(\mathbf{x}, \mathbf{W}(\theta))]$ L-Lipschitz continuous. If there is 243 pair of constants $a > 0$, $b = 0$ such that $\langle \nabla \Phi_F, \nabla \Phi_G \rangle \ge a ||\nabla \Phi_G||_2^2 - b$ and $\varepsilon_t$ meeting 244 requirement 4.2, then the sequence of iterates $\{\mathbf{x}_i\}_{i \in \mathbb{N}}$ generated through the relation

$$246 \quad (4.8) \quad \mathbf{x}_{t+1} = \mathbf{x}_t - \varepsilon_t \nabla_x^C\{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t)$$

247 converges almost surely to a Pareto stationary point.
Proof. Let $B_t = \mathbb{E}[\nabla\Phi_F(\mathbf{x}_t) - \nabla\Phi_G(\mathbf{x}_t)]$. We start from L-Lipschitz smoothness.

$$ \Phi_G(\mathbf{x}_{t+1}) \leq \Phi_G(\mathbf{x}_t) + \langle \nabla\Phi_G(\mathbf{x}_t), \mathbf{x}_{t+1} - \mathbf{x}_t \rangle + \frac{L}{2} \|\mathbf{x}_{t+1} - \mathbf{x}_t\|_2^2. $$

Taking the conditional expectation, and using the definition of SMGD iterations, we have the expression:

$$ \Phi_G(\mathbf{x}_{t+1}) \leq \Phi_G(\mathbf{x}_t) - \mathbb{E}[\varepsilon_t \langle \nabla\Phi_G(\mathbf{x}_t), \nabla_x^C\{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t) \rangle] + \mathbb{E}\left[-\frac{L\varepsilon_t^2}{2} \|\nabla_x^C\{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t)\|_2^2\right] $$

Adding and subtracting $\varepsilon_t \|\Phi_G(\mathbf{x}_t)\|_2^2$, substituting our expression for the bias, and using assumptions 4.2 and 4.3, we see that

$$ \begin{aligned} \Phi_G(\mathbf{x}_{t+1}) &\leq \Phi_G(\mathbf{x}_t) - \varepsilon_t \|\nabla\Phi_G(\mathbf{x}_t)\|_2^2 - \varepsilon_t \langle B_t, \nabla\Phi_G(\mathbf{x}_t) \rangle + \frac{L\varepsilon_t^2}{2}(b + (a_g + 1 - 2a)\|\nabla\Phi_G(\mathbf{x}_t)\|_2^2) \\ &\leq \Phi_G(\mathbf{x}_t) - \varepsilon_t(a - \varepsilon_t \frac{L}{2}M_V)\|\nabla\Phi_G(\mathbf{x}_t)\|_2^2 + \varepsilon_t b + \frac{L\varepsilon_t^2}{2}M_0. \end{aligned} $$

Requiring that $\varepsilon_t < \frac{2a}{LM_V}$ $\forall t$, taking the full expectation, summing, and re-arranging, we get

$$ \sum_{t}^{T} \varepsilon_t \|\nabla \Phi_G(\mathbf{x}_t)\|_2^2 \leq \Phi_G(\mathbf{x}_0) - \mathbb{E}[\Phi_G(\mathbf{x}_T)] + b \sum_{t} \varepsilon_t + \frac{LM_0}{2} \sum_{t} \varepsilon_t^2. $$

Finally, defining $E_T = \sum_t \varepsilon_t$, we can treat $\varepsilon_t/E_T$ as a probability measure. Choosing index $t$ in 1, ..., $T$ and applying Markov's inequality yields

$$ \mathbb{P}(\|\nabla\Phi_G(\mathbf{x})\|_2^2 \geq \epsilon) \leq \frac{\mathbb{E}_{E_T}(\|\nabla\Phi_G(\mathbf{x})\|_2^2)}{\epsilon} \leq \frac{\Phi_G(\mathbf{x}_0) - \mathbb{E}[\Phi_G(\mathbf{x}_T)]}{\epsilon E_T} + \frac{b}{\epsilon} + \frac{LM_0}{2\epsilon E_T} \sum_t \varepsilon_t^2. $$

Taking the limit as $T \to \infty$ gives the result

$$ \mathbb{P}(\|\nabla\Phi_G(\mathbf{x})\|_2^2 \geq \epsilon) \to \frac{b}{\epsilon}. $$

Therefore, if we can find a pair $a, b$ such that $a > 0$ and $b = 0$ in assumption 4.2 then SMGD converges almost surely to a Pareto stationary point. □
Then, even in the case where the bias does not disappear, the SMGD algorithm can be shown to converge to a point on the Pareto set. However, the bias still has an effect. It is straightforward to show that if there is a point, $\mathbf{x}$ where $\mathbb{E}[\nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta))] = \Phi_G(\mathbf{0}) = \mathbf{0}$ and $b = 0$ then SMGD converges to that point.

THEOREM 4.5. Let $F$, $G$, $\varepsilon_t$, $a$, and $b$ be defined as in Thm. 4.4. If there is an $a > 0$ and $b = 0$ meeting the requirements of assumption 4.2 and at least one point, $\mathbf{x}_* \in \mathcal{X}^* := \{\mathbf{x} | \nabla\Phi_F(\mathbf{x}) = \nabla\Phi_G(\mathbf{x}) = \mathbf{0}\}$, then SMGD converges to a point in $\mathcal{X}^*$.

Let $d_t = \mathbf{x}_t - \mathbf{x}^*$ for some $\mathbf{x}^* \in \mathcal{X}^*$ and $B_t = \nabla\Phi_F(\mathbf{x}_t) - \nabla\Phi_G(\mathbf{x}_t)$. To prove Thm. 4.5 it is enough to show that $\langle B_t, d_t \rangle$ goes to zero. Since $B_t \perp d_t$ on the Pareto set implies that $\mathcal{X}^* = \emptyset$ there are two remaining cases. Either $B_t = 0$, indicating that $\nabla\Phi_F(\mathbf{x}_t) = \nabla\Phi_G(\mathbf{x}_t) \to 0$ and that $\mathbf{x}_t \in \mathcal{X}^*$, or $d_t = 0$ in which case $\mathbf{x}_t = \mathbf{x}_*$.

Proof. We start from the inequality

$$d_{t+1}^2 \le d_t^2 - 2\varepsilon_t \langle \nabla_x^C \{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t), d_t \rangle + \varepsilon_t^2 \|\nabla_x^C \{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t)\|_2^2$$

Then, using our definition of $B_t$, we have the expression

$$d_{t+1}^2 \le d_t^2 - 2\varepsilon_t \langle B_t, d_t \rangle + \varepsilon_t^2 (M_v + 1) \|\nabla\Phi_G(\mathbf{x}_t)\|_2^2 + \varepsilon_t^2 M_0.$$

Rearranging and summing we see that

$$(4.9) \qquad \sum_t^T 2\varepsilon_t \langle B_t, d_t \rangle \le d_0^2 - d_T^2 + (M_v + 1) \sum_t^T \varepsilon_t^2 \|\nabla\Phi_G(\mathbf{x}_t)\|_2^2 + M_0 \sum_t^T \varepsilon_t^2.$$

As before, defining $E_T = \sum_t^T \varepsilon_t$, we can use Markov's inequality to show

$$(4.10) \qquad \mathbb{P}_T(\langle B_t, d_t \rangle \ge \epsilon) \le \frac{d_0^2 - d_T^2}{2\epsilon E_T} + (M_v + 1) \frac{\sum_t^T \varepsilon_t^2 \|\nabla\Phi_G(\mathbf{x}_t)\|_2^2}{2\epsilon E_T} + M_0 \frac{\sum_t^T \varepsilon_t^2}{2\epsilon E_T}.$$

Since $\|\nabla\Phi[G](\mathbf{x})\|_2^2 \to 0$ almost surely, we can conclude that the left hand side of Eq. 4.10 goes to zero. □
<figure><img src="image_4.png" /><figcaption>Fig. 3: Comparison of Potential Functions with $f_1(\mathbf{x}, \mathbf{W}(\theta)) = (x - 2 + \mathbf{W}(\theta))^2$ and $f_2(\mathbf{x}, \mathbf{W}(\theta)) = (x + \mathbf{W}(\theta))^2$. SMGDA converges to a point where $\mathbb{E}[\nabla\Phi_F] - \nabla\Phi_G = 0$. Deterministic gradient calculated as $\nabla g_i(\mathbf{x}) = \mathbb{E}[f_i(\mathbf{x}, \mathbf{W}(\theta))].$</figcaption></figure>

**4.2. The Bias of the Stochastic Multigradient.** Why is the stochastic multi-gradient biased when the stochastic gradients used to compute it are not? By assumption we have access to samples of the jacobian of $F$, $J_F(\mathbf{x}, \mathbf{W}(\theta))$ which is an unbiased estimator of the Jacobian of $G$, $J_G(\mathbf{x})$, e.g. $\mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))] = J_G(\mathbf{x})$.

We look to the minimizer $\alpha^*(\mathbf{x}, \mathbf{W}(\theta)) = \operatorname*{argmin}_{\alpha \in \Delta^{k-1}} \alpha^\top J_F(\mathbf{x}, \mathbf{W}(\theta)) J_F(\mathbf{x}, \mathbf{W}(\theta))^\top \alpha$ and notice that it is a quadratic function of the Jacobian. By Jensen's inequality $\mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta)) J_F(\mathbf{x}, \mathbf{W}(\theta))^\top] \ge \mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))] \mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))]^\top$, and so $\mathbb{E}[\alpha^*(\mathbf{x}, \mathbf{W}(\theta))] \ne \alpha^*(\mathbf{x})$.

$\alpha^*$ could be debiased using Monte-Carlo estimation, resampling and averaging the Jacobian and using $\hat{J}_G(\mathbf{x}) = \frac{1}{N} \sum_{i=1}^N J_F(\mathbf{x}, \mathbf{W}(\theta)_i)$ in place of $J_G(\mathbf{x})$ in Eq. 2.7. However, not only would it be computationally taxing to compute $N$ gradients at each iteration, but there would still be residual variance in our estimation rendering our computational efforts somewhat phyrric. Other approaches attempt to smooth the minimizers $\alpha^*$ using exponential smoothing [11]. This approach effectively reduces the variance of $\alpha^*$ and enforces that the estimator $\alpha^*$ be smooth across iterations. In contrast, our approach to the problem is to mitigate the bias at its source, in the calculation of $\alpha^*$ itself.

**5. Proposed Approach: Stochastic Tangential Pareto Dynamics.** The SMGD algorithm, even if run from distinct starting points, converges to a subset of the Pareto set. Then, even if completely debiased, in order to estimate the whole Pareto
set one has to restart the algorithm from several points in design space. Not only is this approach not guaranteed to uncover the whole Pareto set but it is computationally taxing. To ameliorate this, we propose to modify the stochastic multigradient descent iterations in two ways. First, we introduce a novel debiasing strategy for the stochastic multigradient. This allows the algorithm to converge to the whole of the Pareto front. Then, we add a noise term which guides the dynamics generated by the algorithm tangentially along the Pareto set, promoting exploration as well as overcoming any residual bias. Because of the added noise term, this approach efficiently explores the area around the Pareto set and is less likely to stay stuck in a saddle point. Furthermore, iterations can be performed with single samples of the quantities of interest and their gradients at each round.

Dubbed stochastic tangential Pareto dynamics (STPD), the points generated by this approach create a dense sample set concentrated on the Pareto set.

**5.1. Debiasing the Stochastic Multigradient.** We would like to use the deterministic Jacobian matrix, $J_G$, to calculate the multigradient. This would be equivalent to having direct access to $\mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))]\mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))]^\top$. However, we only have access to samples of the matrix $J_F(\mathbf{x}, \mathbf{W}(\theta))J_F(\mathbf{x}, \mathbf{W}(\theta))^\top$. We notice that the two quantities are related,

(5.1)

$$ \mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))J_F(\mathbf{x}, \mathbf{W}(\theta))^\top] = \Sigma_{J_F(\mathbf{x}, \mathbf{W}(\theta))} + \mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))]\mathbb{E}[J_F(\mathbf{x}, \mathbf{W}(\theta))]^\top, $$

where $\Sigma_{J_F(\mathbf{x}, \mathbf{W}(\theta))}$ denotes the covariance matrix of the gradients of $F$.

We propose to debias the stochastic multigradient using exponential smoothing. We first compute an online estimate of $\Sigma_{J_F(\mathbf{x}, \mathbf{W}(\theta))}$ using the recursive estimate

$$ (5.2) \qquad \hat{\mu}_{t+1} = \gamma_t \hat{\mu}_t + (1 - \gamma_t) J_F(\mathbf{x}_t, \mathbf{W}(\theta)_t) $$

$$ (5.3) \qquad \hat{\Sigma}_{J_F, t+1} = \gamma_t \hat{\Sigma}_{J_F, t} + (1 - \gamma_t) (J_F(\mathbf{x}_t, \mathbf{W}(\theta)_t) J_F(\mathbf{x}_t, \mathbf{W}(\theta)_t)^\top - \hat{\mu}_t \hat{\mu}_t^\top) $$

using individual samples of $J_F(\mathbf{x}, \mathbf{W}(\theta))$ and $\gamma_t \in (0, 1)$. We then calculate $\alpha^*$ using
the following modified formulation.

$$ (5.4) \quad \hat{\alpha}^*(\mathbf{x}) = \underset{\alpha \in \Delta^{k-1}}{\operatorname{argmin}} \alpha^\top (J_F(\mathbf{x}, \mathbf{W}(\theta)) J_F(\mathbf{x}, \mathbf{W}(\theta)))^\top - \hat{\Sigma}_{J_F(\mathbf{x}, \mathbf{W}(\theta))} \alpha. $$

We can then calculate a debiased form of the stochastic gradient

$$ (5.5) \quad \widehat{\nabla_x^C \{F\}}(\mathbf{x}, \mathbf{W}(\theta)) := \sum_i \hat{\alpha}_i^* \nabla f_i(\mathbf{x}, \mathbf{W}(\theta)) $$

This approach has the advantage of not requiring excess computation in individual rounds and places less assumptions on the required behavior of $\alpha^*$ while also effectively reducing the bias in the calculation of $\nabla_x^C \{F\}(\mathbf{x}, \mathbf{W}(\theta))$.

**5.2. Stochastic Tangential Pareto Dynamics.** To find and explore the Pareto set we generate a sequence of points using the recurrence relation,

$$ (5.6) \quad \mathbf{x}_{t+1} = \mathbf{x}_t - \varepsilon_t \widehat{\nabla_x^C \{F\}}(\mathbf{x}_t) + \sqrt{2\varepsilon_t \beta_t} \Pi_t Z, $$

with $\nabla_x^C \{F\}(\mathbf{x}_t)$ a direction of descent, $Z \sim \mathcal{N}(0, \mathbf{1}_{d \times d})$, $\Pi_t$ a projection operator that projects the noise perpendicular to the direction of descent, $\varepsilon_t$ a stepsize parameter, and $\beta_t$ the scale of the noise. The subscript $t$ in this case should be understood as being adapted to the filtration generated by $\mathbf{x}_t$. This algorithm can be intuitively understood as gradient descent with an additive exploratory noise term. To keep the samples from converging to a limiting distribution centered on a single point we set the stepsize parameter, $\varepsilon_t$, and the moving average parameter, $\gamma_t$ as a pair of nonincreasing but strictly positive numbers:

$$ (5.7) \quad \begin{aligned} \{\varepsilon_t\}_{t=0,...,\infty} &:= \varepsilon_0 \ge \varepsilon_1 \ge \dots \ge \varepsilon_\infty > 0, \\ \{\gamma_t\}_{t=0,...,\infty} &:= \gamma_0 \ge \gamma_1 \ge \dots \ge \gamma_\infty > 0. \end{aligned} $$

The moving average parameter, as we see below, is used in the estimation of the debiased stochastic multigradient and its expected value.

Ideally, we would have access to true values of the gradients for each objective.
It would then be possible to use the deterministic multigradient, $\nabla_x^C\{G\}(\mathbf{x})$, as an unbiased direction of descent. However, having only samples of the gradients of each objective, we approximate the direction of descent. Here, the expected value of the debaised stochastic multigradient as given in Eq. 5.5 defines the gradient of the potential, $\nabla\Phi_F(\mathbf{x})$. Since we do not have access to its expected value, we approximate it using single samples.

The performance of the estimation of the debiased stochastic multigradient depends on the noise in the stochastic gradients and the strength of the additive noise term. If the noise in the individual stochastic gradients is low then exponential smoothing with a large value of $\gamma_\infty$ will yield a good approximation of the gradient. As either the strength of the additive noise term or the noise in the gradients increases the autocorrelation between iterates falls and the exponential smoothing procedure breaks down. This is viewed with more detail in our companion paper.

$\Pi_t$ removes the component of the additive noise that lies parallel to the multigradient. On the Pareto set its value is unity and the algorithm produces iterates that move randomly. Near the Pareto set, the direction orthogonal to the multigradient lies tangential to the Pareto set. The tangential motion of the iterates generated by STPD can then be understood as a tangential drift along the Pareto set. If $\Pi_t$ was set permanently to unity then this approach would degrade into a type of stochastic Langevin dynamics. This would be unsatisfactory, however, as the variance of samples orthogonal to the Pareto set would be much higher.

To calculate $\Pi_t$ we need access to the expected value of the direction of descent. However, we do not in general have access to an analytic expression for the expected value of the debiased multigradient. Therefore it must be estimated during the course of optimization. This estimation is also biased due to the nonzero stepsize.

At each iteration we update the expected value of the stochastic multigradient

$$ (5.8) \qquad \overline{\nabla_x^C\{F\}}_{t+1} = \gamma_t \overline{\nabla_x^C\{F\}}_t + (1-\gamma_t) \overline{\nabla_x^C\{F\}}(\mathbf{x}_t). $$

To calculate the projection operator we use the exponentially smoothed debiased
stochastic multigradient,

$$ (5.9) \qquad \Pi_{t,i,j} = \delta_{i,j} - \frac{1}{\|\overline{\nabla_x^C \{F\}}_t\|^2} \overline{\nabla_x^C \{F\}}_{t,i} \overline{\nabla_x^C \{F\}}_{t,j}, $$

where $\delta_{i,j}$ denotes the dirac delta function.

The value of $\beta_t$ determines the relative strength of the additive noise compared to the direction of descent and can be viewed as a type of temperature parameter. We only require that $\beta_t$ be a positive upper bounded function adapted to the filtration generated by $\mathbf{x}_t$. It must be positive since a value of zero would remove the noise term, and it must be upper bounded in order to allow for convergence. We discuss possible function choices in Section 5.4. The effect of the relative strength of the noise term is studied numerically in our companion work.

This entire process is summarized in algorithm 5.1.

<caption>**Algorithm 5.1 Summary of Stochastic Tangential Pareto Dynamics**</caption>

Input {$\varepsilon_t$}, {$\gamma_t$}, {$\mathbf{x}_0$}.
Initialize $\hat{\mu}_0 = \mathbf{0}$, $\hat{\Sigma}_{J_F,0} = \mathbf{1}_{d \times d}$, $t = 1$, $\overline{\nabla_x^C \{F\}}_0 = 0$.
**while** Running **do**
Query Stochastic Oracle for $J_F(\mathbf{x}_t, \mathbf{W}(\theta)_t)$ ▷ Referred to as $J_{F_t}$ below.
$\hat{\mu}_t \leftarrow \gamma_t \hat{\mu}_{t-1} + (1 - \gamma_t) J_{F,t}$
$\hat{\Sigma}_{J_F,t} \leftarrow \gamma_t \hat{\Sigma}_{J_F,t-1} + (1 - \gamma_t)(J_{F,t} - \hat{\mu}_t)(J_{F,t} - \hat{\mu}_t)^\top$
$\hat{\alpha}^* \leftarrow \underset{\alpha \in \Delta^{k-1}}{\operatorname{argmin}} \alpha^\top (J_{F,t} J_{F,t}^\top - \hat{\Sigma}_{J_F,t}) \alpha$ ▷ Eq. ??
$\overline{\nabla_x^C \{F\}}_t \leftarrow J_{F,t}^\top \hat{\alpha}^*$ ▷ Eq. ??
$\overline{\nabla_x^C \{F\}}_t \leftarrow \gamma_t \overline{\nabla_x^C \{F\}}_{t-1} + (1 - \gamma_t) \overline{\nabla_x^C \{F\}}_t$ ▷ Eq. 5.8
$\Pi_t \leftarrow \mathbf{1}_{d \times d} - \frac{\overline{\nabla_x^C \{F\}}_t \overline{\nabla_x^C \{F\}}_t^\top}{\|\overline{\nabla_x^C \{F\}}_t\|^2} \overline{\nabla_x^C \{F\}}_t \overline{\nabla_x^C \{F\}}_t^\top$ ▷ Eq. 5.9
$Z \sim \mathcal{N}(\mathbf{0}, \mathbf{1}_{d \times d})$
$\mathbf{x}_{t+1} \leftarrow \mathbf{x}_t - \varepsilon_t \overline{\nabla_x^C \{F\}}_t + \sqrt{2\varepsilon_t \beta_t} \Pi_t Z$ ▷ Eq. 5.6
$t \leftarrow t + 1$
**end while**

396

5.3. **Stochastic Tangential Pareto Dynamics Samples the Whole Pareto**
**Front.** Here we show that, given enough time, stochastic tangential Pareto dynamics samples arbitrarily close to all points on the Pareto set.

In addition to assumption 4.2 and L-Lipschitz smoothness of $\Phi_G$, we assume coerciveness of $\Phi_G$.

401
402 ASSUMPTION 5.1. ∃c_a, c_b ∈ ℝ+, ∀x, such that

$$ (5.10) \qquad ||\nabla \Phi_G(\mathbf{x})||^2 \ge c_a \Phi_G(\mathbf{x}) - c_b, \qquad ||\mathbf{x}||^2 \le c_a \Phi_G(\mathbf{x}) + c_b . \qquad 403 \qquad 404 $$

This assumption is coercive in the sense that the strength the gradient increases the further away $\mathbf{x}$ is from the minima, "pushing" iterates back towards the minima as the value of the scalar potential gets larger. Assumption 5.1 is necessary to show the ergodicity of the SDE $dX = \nabla F(X) + dW$, and so it is appropriate in this context as well [10, 8, 2].

410 We can now show the main theorem.

411 THEOREM 5.2. Given a set of constants {$\epsilon_t$}$_{t \in \mathbb{N}}$, assuming $\Phi_G : \mathbb{R}^d \to \mathbb{R}$ is L-Lipschitz smooth and obeys conditions 5.1 and 4.2 almost everywhere, and assuming lower bound $B := L\beta^{-1}d + b + \frac{\epsilon_0 LM_0}{2} + d + c_b \le c_1 \inf \Phi_G(\mathbf{x})$, then for any $\mathbf{y}^* \in P^{-1}$ and $\epsilon, p > 0$

$$ (5.11) \qquad \mathbb{P}(||\mathbf{x}_t - \mathbf{y}^*||| \le \epsilon \text{ for some } t < \infty) \ge 1 - p. \qquad 415 $$

Our framework for proving Theorem 5.2 for STPD follows closely from [8, 2, 10]. We first show recurrence, that there is a compact sublevel set, $\Phi_G(\mathbf{x}) \le M$, that is reached infinitely many times by STPD. We then show reachability, on a compact sublevel set, there is a nonzero probability to reach any target point $\mathbf{y}^* \in \mathcal{P}$. Once we have established both recurrence and reachability, it is straightforward to see that STPD reaches an arbitrary point $\mathbf{y}^*$ in finite time with probability arbitrarily close to 1. Since $\mathbf{y}^*$ is arbitrary, it follows that it reaches the whole Pareto set.

423 First we show recurrence.

424 LEMMA 5.3. Recurrent visits to the sublevel set M.

425 Assume that $\Phi_G$ is L-Lipschitz smooth, obeys assumptions 5.1 and 4.2, and let {$\epsilon_t$}$_{t \in \mathbb{N}}$ be a nonincreasing sequence such that $\epsilon_0 \ge \epsilon_1, ..., \ge \epsilon_\infty > 0$. Let $M > 0$ and define the constant $B := L\beta^{-1}d + b + \frac{\epsilon_0 LM_0}{2} + c_b \le c_a \inf \Phi_G(\mathbf{x})$. Given a sequence of stopping
428 times $\tau_{k+1} = \inf\{t: t > \tau_k, \Phi_G(\mathbf{x}_t) \le M\}$ then

$$429 \quad (5.12) \qquad \text{a)} \quad \tau_0 = \frac{\log\left[\left[\frac{\Phi_G(\mathbf{x}_0)}{M - \frac{\varepsilon_0 B}{\varepsilon_\infty a}}\right]^2\right]}{2c_a \varepsilon_\infty}$$

430 and

$$431 \quad (5.13) \qquad \text{b)} \quad \mathbb{E}[\tau_j] = \tau_0 + (j+1)M$$

432 For proof see the appendix.

433 We then have reachability.

434 LEMMA 5.4 (Reachability). Given that $\Phi_G(\mathbf{x}_{\tau_k}) \le M$, let $\mathbb{E}[\nabla\Phi_F(\mathbf{x})] \le \mathbf{D}$ and
435 $c_t > 0$. We have, for $\mathbf{y}^* \in \{\mathbf{y} : \Phi_G(\mathbf{y}) \le M\} \cap \{\mathbf{y} : \mathbf{y} \in \mathbf{P}^{-1}\}$

$$436 \quad (5.14) \qquad \mathbb{P}(|\mathbf{x}_{\tau_j+t} - \mathbf{y}^*|^2 \le \epsilon) = c_t > 0$$

437 For proof see the appendix.

438 Having shown both recurrance and reachability established in lemmas 5.3 and
439 5.4, we now prove Theorem 5.2

440 Proof. First, define a sequence of stopping times and a $\tau_{j+1} = \inf\{t: t >$
441 $\tau_j, \Phi_G(\mathbf{x}_{\tau_j+1}) \le M\}$. Then, also define a stopping time $\tau^* = \inf t: ||\mathbf{x}_t - \mathbf{y}^*|| \le \epsilon$
442 for $\epsilon > 0$. We have that

$$443 \quad \mathbb{P}(\tau^* \ge T) = \mathbb{P}(\tau^* \ge T, \tau_j \ge T) + \mathbb{P}(\tau^* \ge T, \tau_j \le T) \\
444 \quad = \mathbb{P}(\tau_j \ge T) + \mathbb{P}(\tau^* \ge T, ||\mathbf{x}_{\tau_k} - \mathbf{y}^*|| \ge \epsilon \quad \forall k = \{1, \dots, j\}) \\
445 \quad \le \mathbb{E}[\tau_j] + \prod_{k=1}^j \mathbb{P}(|\mathbf{x}_{\tau_k} - \mathbf{y}^*| \ge \epsilon) \quad (\text{Lem. 5.3})$$

$$446 \quad \le \frac{(j+1)M + \tau_0}{T} + \prod_{k=1}^{j} (1 - c_j) \quad (\text{Lem. 5.4})$$

$$447 \quad \le \frac{(j+1)M + \tau_0}{T} + (1 - c^*)^j. \\ 448$$
Where we have used the boundedness of $c_j$ in the final line. Setting $\frac{(j+1)M+\tau_0}{T} = \alpha p$ and $(1-c^*)^j = (1-\alpha)p$ for $\alpha \in (0, 1)$ and solving for $j$ and $T$ we see that setting

$$j \ge \frac{\ln((1-\alpha)p)}{1-c^*} \quad (\text{and}) \quad T \ge \frac{\frac{\ln((1-\alpha)p)}{1-c^*} + 1}{\alpha p} M + \tau_0$$

gives us $\mathbb{P}(\tau^* \ge T) \le p$. Since this is true for any $p$,

$$\mathbb{P}(\tau^* \le T) \ge 1 - p,$$

which can be made arbitrarily close to 1.

**5.4. Discussion.** Here we consider implementation choices, the effect of non-stationarity of the sample distribution, parameter settings, and dealing with saddle points and nonconvexity.

We have presented here a basic version of our approach with the aim of imparting understanding about the general working of the algorithm. There are other techniques which may be included in the algorithm. They do not affect our theoretical results and may improve the sampling properties of the algorithm in practice.

We first consider variance reduction. Variance reduction techniques can be used on both $\nabla f_i(\mathbf{x}, \mathbf{W}(\theta))$ and the debiased stochastic multigradient $\nabla_x^C\{\mathbf{F}\}(\mathbf{x})$. Reducing the variance of the objective gradients would further decrease the bias of the calculation of $\nabla_x^C\{\mathbf{F}\}(\mathbf{x})$. It would also cause the points generated by the algorithm to be more concentrated on the Pareto set. Common approaches include the popular SAGA, SVRG, and minibatching [1]. The increased precision comes at a cost, for a given stepsize lower gradient variance will also decrease the speed at which the algorithm diffuses along the Pareto set. The use of variance reduction is ultimately problem dependent and will affect parameter choices for stepsize and additive noise strength.

It is also possible to include preconditioners in our approach, i.e. take steps

$$\mathbf{x}_{t+1} = \mathbf{x}_t - \varepsilon_t H_t \Phi_F(\mathbf{x}) + \sqrt{2\beta_t \varepsilon_t H_t} \Pi_t Z$$
where $H_t$ is a positive definite matrix. There are many options for the preconditioner, as discussed in the context of single objective stochastic Langevin dynamics in [7].

The samples generated by STPD do not reach a stationary distribution. This is because of the nonvanishing stepsize. As a result the estimates for the jacobian, its covariance, and the debiased stochastic multigradient do not converge. This appears in the algorithm as a bias in both the direction of descent and the projection operator. This bias is ultimately overcome by the variation introduced by the additive noise.

The bias in the direction of descent lies tangential to the Pareto set. Clearly, immediately projecting normal to this direction would be counterproductive. Averaging over the history of the direction of descent allows us to estimate a true tangential direction. This average, however, still contains a trace of the bias, causing the projected noise to have a component normal to the Pareto set. This manifests as an increased residual variance of the samples around the Preto set.

Projecting normal to the average history of the direction of descent is not the only way to calculate the projection operator. It is also possible to average the history of the directions perpendicular to the direction of descent, e.g. calculate $\Pi_t$ as:

$$ \Pi_{t+1} = \gamma_t \Pi_t + \left(1 - \frac{1}{\|\widehat{\nabla_x^C \{F\}}(\mathbf{x}_t)\|_2^2} \widehat{\nabla_x^C \{F\}}(\mathbf{x}_t) \widehat{\nabla_x^C \{F\}}(\mathbf{x}_t)^\top\right) $$

We have chosen exponential smoothing to model several quantities. It has the advantage of requiring little storage and being easy to calculate. However, its effectiveness depends strongly on the correlation time of the sample gradients. The correlation time, in turn, is influenced by the noise of the stochastic gradients along the Pareto set and the relative strength of the additive noise term. These interactions are not straightforward. More complex models can be used in the place of simple exponential smoothing. However, they will similarly show a dependence on the correlation time of the sample gradients and must be calibrated. This can be done simply by minimizing the predictive error of the estimation of a quantity of interest, e.g. the Jacobian.

It has been previously mentioned that the additive noise term $\beta_t$ need only be a
positive upper bounded function adapted to the filtration generated by $\mathbf{x}_t$. This leaves
flexibility in its choice beyond a constant function. Since the additive noise term may
interfere with initial convergence of the iterates to the Pareto set it might be beneficial
to damp it when iterates are far from the Pareto set. The strength of the additive noise
can be damped far from the Pareto set by selecting $\beta_t = C e^{-\overline{\nabla_x^c \{F\}}^\top} A \overline{\nabla_x^c \{F\}} + D$ with
$C, D \in \mathbb{R}^{d \times d}$, $A \in \mathbb{R}_+^{d \times d}$. The dynamics generated by STPD are also strongly affected
by the relative strength of the additive noise term. This suggests a scaling to augment
the strength of the additive noise in proportion to the variation of the objective
functions by selecting $\beta_t = \sum_i \text{VAR}(\nabla f_i(\mathbf{x}_t, \mathbf{W}(\theta)_t))$. Vector valued functions can
also be considered.

We take a moment to discuss saddle points in the scalar potential. We have
proven that the iterates generated by STPD cover the whole Pareto set under mild
assumptions. This does not preclude cases with non-convex objectives. It has been
shown in the single objective case that adding Gaussian noise to stochastic gradients
can allow gradient descent algorithms to escape shallow second order minima [2]. In
the multi-objective case, local minima in the objective functions correspond to saddle
points in the scalar potential which, if mild, may be similarly escaped. However,
a potential with a split Pareto set (a 'W' shaped potential) may lead to complica-
tions. Either multiple restarts or an alternative approach more suited to nonconvex
optimization may be more effective for such situations.

Diffusing across the Pareto set is not the only way to exploit the stochastic multi-
gradient. An alternative approach would be to select an ensemble of random points
in design space and optimize them individually using debiased SMGDA. This ap-
proach may be more effective on non-convex problems and problems with discontin-
uous Pareto sets. However, it has two corresponding drawbacks. One is the loss of
information. Two points from distinct regions in design space may converge to similar
locations on the Pareto set, making one redundant. Another drawback is resolution.
Picking an ensemble of points pre-determines the resolution at which one can deter-
mine the Pareto set. This has a similar flavor to picking the number of iterations of
STPD. The computational cost, however, is not the same. STPD scales linearly with
the number of samples taken, whereas adding another point in an ensemble algorithm scales multiplicatively with the number of iterations used.

**6. Conclusion.** We have introduced a new approach to solving multi-objective optimization problems using noise added unbiased stochastic multigradients, in which the whole of the Pareto front is of interest, stochastic tangential Pareto dynamics. Along the way, we have reformulated our multi-objective optimization problem as the minimization of an equivalent scalar potential function. Using this reformulation, we have presented alternative proofs of convergence for the MGD algorithm along with new proofs of convergence for the SMGD algorithm; showing that, despite bias, the SMGD algorithm converges to a Pareto stationary point in the convex, strongly convex, and nonconvex L-Lipshitz smooth case. We have also shown that, when there are points on the Pareto set where the stochastic multigradient is unbiased, that SMGDA converges to one of those points. Finally, we have shown that STPD provably generates samples arbitrarily close to the whole of the Pareto set. Our approach produces a noisy snapshot of the Pareto front and reducing this variance through alternative gradient estimators and/or second order algorithms remains a promising avenue of future research.

A complementary line of inquiry involves characterizing the Pareto front, and set, from a fixed set of samples. In our companion paper we give a probabilistic characterization of the Pareto set across the whole of design space, giving a straightforward recipe to not only infer the Pareto set from STPD but to postprocess the results of any stochastic sampling algorithm that yields samples concentrated along the Pareto set [CITATION NEEDED]. Taken in conjunction with STPD this yields a straightforward way to estimate all likely members of the Pareto set of a stochastic program.

**Acknowledgments.** We would like to acknowledge ourselves for being spectacular.
561 [1] L. BOTTOU, F. E. CURTIS, AND J. NOCEDAL, Optimization methods for large-scale machine
learning, 2018, https://arxiv.org/abs/1606.04838.

563 [2] X. CHEN, S. S. DU, AND X. T. TONG, On stationary-point hitting time and ergodicity of
stochastic gradient langevin dynamics, 2020, https://arxiv.org/abs/1904.13016.

565 [3] M. FRANK AND P. WOLFE, An algorithm for quadratic programming, Naval Research Logistics
Quarterly, 3 (1956), pp. 95–110, https://api.semanticscholar.org/CorpusID:122654717.

567 [4] Z. HU, K. SHALOUDEGI, G. ZHANG, AND Y. YU, Federated learning meets multi-objective opti-
mization, 2023, https://arxiv.org/abs/2006.11489, https://arxiv.org/abs/2006.11489.

569 [5] S. LIU AND L. N. VICENTE, The stochastic multi-gradient algorithm for multi-objective opti-
mization and its application to supervised machine learning, 2021, https://arxiv.org/abs/
1907.04472.

572 [6] M. M. MÄKELÄ, V.-P. ERONEN, AND N. KARMITSA, On Nonsmooth Multiobjective Opti-
mality Conditions with Generalized Convexities, Springer New York, New York, NY,
2014, pp. 333–357, https://doi.org/10.1007/978-1-4939-0808-0_17, https://doi.org/10.
1007/978-1-4939-0808-0_17.

576 [7] G. MARCEAU-CARON AND Y. OLLIVIER, Natural langevin dynamics for neural networks, in
Geometric Science of Information, F. Nielsen and F. Barbaresco, eds., Cham, 2017, Springer
International Publishing, pp. 451–459.

579 [8] J. MATTINGLY, A. STUART, AND D. HIGHAM, Ergodicity for sdes and approximations: locally
lipschitz vector fields and degenerate noise, Stochastic Processes and their Applications,
101 (2002), pp. 185–232, https://doi.org/https://doi.org/10.1016/S0304-4149(02)00150-3,
https://www.sciencedirect.com/science/article/pii/S0304414902001503.

583 [9] Q. MERCIER, Optimisation multicritère sous incertitudes : un algorithme de descente stochas-
tique, theses, COMUE Université Côte d'Azur (2015 - 2019), Oct. 2018, https://theses.hal.
science/tel-02063322.

586 [10] S. MEYN AND R. TWEEDIE, Markov Chains and Stochastic Stability, vol. 92, 01 1993, https:
//doi.org/10.2307/2965732.

588 [11] S. ZHOU, W. ZHANG, J. JIANG, W. ZHONG, J. GU, AND W. ZHU, On the con-
vergence of stochastic multi-objective gradient manipulation and beyond, in Ad-
vances in Neural Information Processing Systems, S. Koyejo, S. Mohamed, A. Agar-
wal, D. Belgrave, K. Cho, and A. Oh, eds., vol. 35, Curran Associates,
Inc., 2022, pp. 38103–38115, https://proceedings.neurips.cc/paper_files/paper/2022/file/
f91bd64a3620aad8e70a27ad9cb3ca57-Paper-Conference.pdf.

**Appendix A. Proofs of Recurrence and Reachability of STPD.** In order
to show that STPD samples the whole of the Pareto set we must show both recurrence
and reachability. Recurrence shows that the algorithm iterates reach a sublevel set
597 upper bounded by $M$ an infinite number of times. Reachability shows that, having
598 reached the sublevel set upper bounded $M$, there is a nonzero probability of reaching
599 any point on the Pareto set inside of that sublevel set. We begin with the proof of
600 recurrence.

601 LEMMA. Recurrent visits to the sublevel set $M$.

602 Assume that $\Phi_G$ is L-Lipschitz smooth, obeys assumptions 5.1 and 4.2, and let $\{\varepsilon_t\}_{t \in \mathbb{N}}$
603 be a nonincreasing sequence such that $\varepsilon_0 \ge \varepsilon_1, \dots, \ge \varepsilon_\infty > 0$. Let $M > 0$ and define
604 the constant $B := L \max(\beta_t)d + b + \frac{\varepsilon_0 LM_0}{2} + c_b \le c_a \inf \Phi_G(\mathbf{x})$. Given a sequence of
605 stopping times $\tau_{k+1} = \inf\{t : t > \tau_k, \Phi_G(\mathbf{x}_t) \le M\}$ then

$$
606 \quad a) \quad \tau_0 = \frac{\log \left[ \left[ \frac{\Phi_G(\mathbf{x}_0)}{M - \frac{\varepsilon_0 B}{\varepsilon_\infty a}} \right]^2 \right]}{2c_a \varepsilon_\infty}
$$

607 and

$$
608 \quad b) \quad \mathbb{E}[\tau_j] = \tau_0 + (j+1)M
$$

609
Proof of lemma 5.12(a) (Recurrence). Let $\zeta_t = \langle \nabla \Phi_G(\mathbf{x}_t), \nabla \Phi_F(\mathbf{x}_t) - \nabla \Phi_G(\mathbf{x}_t) \rangle$.

610 We start from L-Liptshitz smoothness of $\Phi_G$.

$$
611 \quad \Phi_G(\mathbf{x}_{t+1}) \le \Phi_G(\mathbf{x}_t) + \langle \nabla \Phi_G(\mathbf{x}_t), \mathbf{x}_{t+1} - \mathbf{x}_t \rangle + \frac{L}{2} \|\mathbf{x}_{t+1} - \mathbf{x}_t\|^2. \\
612
$$

613 Substituting in $\nabla_x^c\{F\}(\mathbf{x})$, adding and subtracting $\varepsilon_t ||\Phi_G(\mathbf{x}_t)||^2$ and using 4.7,

$$
614 \quad \le \Phi_G(\mathbf{x}_t) - \varepsilon_t ||\nabla \Phi_G(\mathbf{x}_t)||^2 - \varepsilon_t \zeta_t + \frac{L}{2} (M_{V_0} + M_V ||\nabla \Phi_G(\mathbf{x}_t)||^2) \\
615 \quad \qquad + \sqrt{2\varepsilon_t \beta_t} \langle \nabla \Phi_G(\mathbf{x}_t), \Pi_t Z \rangle + \frac{L}{2} (2\sqrt{2\varepsilon_t \beta_t} \langle \nabla \Phi_F(\mathbf{x}_t), \Pi_t Z \rangle + 2\varepsilon_t \beta_t \Pi_t \Pi_t^\top \langle Z, Z \rangle).
$$

616
617
618 Taking the conditional expectation and simplifying,

$$619 \quad \mathbb{E}_t[\Phi_G(\mathbf{x}_{t+1})|\mathbf{x}_t] \le \Phi_G(\mathbf{x}_t) - \varepsilon_t(a - \frac{\varepsilon_t L(M_v+1)}{2})||\nabla\Phi_G(\mathbf{x}_t)||^2 + \varepsilon_t(b + \frac{L\varepsilon_t}{2}M_0 + L\beta_t d)$$

$$620 \quad \le \Phi_G(\mathbf{x}_t) - \varepsilon_t ||\nabla\Phi_G(\mathbf{x}_t)||^2 + \varepsilon_t \left(b + \frac{L\varepsilon_t}{2}M_0 + L\beta_t d\right)$$

$$621 \quad \le (1 - c_a \varepsilon_t) \Phi_G(\mathbf{x}_t) + \varepsilon_t B$$

$$623 \quad \le e^{-c_a \varepsilon_t} \Phi_G(\mathbf{x}_t) + \varepsilon_t B.$$

624 Where we have used the inequality $1 - x \le e^{-x}$ and enforced the requirement that

625 $\varepsilon_0 \le \frac{2a}{L(M_v+1)}$. Taking the full expectation, iterating, and setting the above less than

626 equal to $M$ gives the relation,

$$627 \quad \mathbb{E}[\Phi_G(\mathbf{x}_{\tau_0})] \le e^{-c_a \tau_0 \varepsilon_\infty} \Phi_G(\mathbf{x}_0) + \frac{\varepsilon_0 B}{\varepsilon_\infty a} \le M.$$

628 Solving for $\tau_0$ gives the result. □

629 For the Proof of part b of the lemma, we will first prove that the quantity

$$630 \quad \Phi_G(\mathbf{x}_{t \wedge \tau_{j+1}}) + t \wedge \tau_{j+1}$$

631 is a supermartingale with respect to $\tau_j$.

LEMMA A.2.

$$632 \quad \Phi_G(\mathbf{x}_{t \wedge \tau_{j+1}}) + t \wedge \tau_{j+1}$$

633 is a supermartingale with respect to $\tau_j$.

634 *Proof.* Since $t \wedge \tau_{j+1}$ is a stopping time and a martingale, it suffices to show that

$$635 \quad \mathbb{E}_{\tau_j} [\Phi_G(\mathbf{x}_{\tau_{j+1}})] \le \Phi_G(\mathbf{x}_{\tau_j})$$
Note that

$$637 \qquad \mathbb{E}_{\tau_j}[\Phi_G(\mathbf{x}_{\tau_{j+1}})] \le (1 - a\varepsilon_{\tau_j})\Phi_G(\mathbf{x}_{\tau_j}) + \varepsilon_{\tau_j}B \le M$$

Is true as we can always pick $a \ge 1$. Since $\Phi_G(\mathbf{x}_{\tau_j}) \ge \inf \Phi_G(()) \ge \frac{B}{a}$,

$$639 \qquad \mathbb{E}_{\tau_j}[\Phi_G(\mathbf{x}_{\tau_{j+1}})] \le \Phi_G(\tau_j)$$

And we have the result.

□

Having proved lemma A.2 we can go on to prove part b of lemma 5.3

Proof of 5.4(b) (Recurrence). Since

$$643 \qquad \mathbb{E}[\Phi_G(\mathbf{x}_{t \wedge \tau_{j+1}})|\tau_j] + t \wedge \tau_{j+1} \le \Phi_G(t \wedge \tau_j) + t \wedge \tau_j$$

we can allow $t \to \infty$ to see that

$$645 \qquad \mathbb{E}[\tau_{j+1}|\tau_j] \le \mathbb{E}[\Phi_G(\mathbf{x}_{\tau_{j+1}})|\tau_j] + \tau_{j+1} \le \Phi_G(\mathbf{x}_{\tau_j}) + \tau_j$$

Iterating, using the relation $\Phi_G(\mathbf{x}_{\tau_j}) \le M$, and taking the full expectation we see that

□

$$648 \qquad \mathbb{E}[\tau_{j+1}] \le (j+1)M + \tau_0$$

We now prove reachability.

REACHABILITY. Given that $\Phi_G(\mathbf{x}_{\tau_k}) \le M$, let $\mathbb{E}[\nabla\Phi_F(\mathbf{x})] \le \mathbf{D}$ be the componentwise maximum value of the gradient in the sublevel set, and $c_t > 0$. We have, for $\mathbf{y}^* \in \{\mathbf{y} : \Phi_G(\mathbf{y}) \le M\} \cap \{\mathbf{y} : \mathbf{y} \in \mathbb{P}^{-1}\}$

$$653 \qquad \mathbb{P}(|\mathbf{x}_{\tau_j+t} - \mathbf{y}^*|^2 \le \epsilon) = c_t > 0$$
Proof of Lemma 5.4 (Reachability). for a finite t we have the events

$$ \mathcal{A} = \|\mathbf{x}_{\tau_j+t} - \mathbf{y}^*\|^2 \le \epsilon \quad \mathcal{B} = \|\mathbf{x}_{\tau_j} - \mathbf{y}^* + \sum_{s=\tau_j}^{\tau_j+t-1} \epsilon_s \mathbf{D} + \sum_{s=\tau_j}^{\tau_j+t-1} \sqrt{2\epsilon_s \beta_s} \Pi_s Z_s\|^2 \le \epsilon $$

with $Z \sim \mathcal{N}(0, \mathbb{1}_{d \times d})$. Since $t < \infty$ we conclude that $\mathbb{P}(\mathcal{B}) > 0$. From the definition of STPD, we have

$$ \begin{aligned} \mathbf{x}_{\tau_j+t} &= \mathbf{x}_{\tau_j} - \sum_{s=\tau_j}^{\tau_j+t-1} \nabla \Phi_F(\mathbf{x}_s) + \sum_{s=\tau_j}^{\tau_j+t-1} \sqrt{2\epsilon_s \beta_s} \Pi_s Z_s \\ &\le \mathbf{x}_{\tau_j} + \sum_{s=\tau_j}^{\tau_j+t-1} D + \sum_{s=\tau_j}^{\tau_j+t-1} \sqrt{2\epsilon_s \beta_s} \Pi_s Z_s \end{aligned} $$

Subtracting $\mathbf{y}^*$ from both sides and taking the norm, we see that

$$ \|\mathbf{x}_{\tau_j+t} - \mathbf{y}^*\|^2 \le \|\mathbf{x}_{\tau_j} - \mathbf{y}^* + \sum_{s=\tau_j}^{\tau_j+t-1} D + \sum_{s=\tau_j}^{\tau_j+t-1} \sqrt{2\epsilon_s \beta_t} \Pi_s Z_s\|^2 $$

Where we can see that on the left hand side we have event $\mathcal{A}$. Since event $\mathcal{A}$ occurs almost surely if event $\mathcal{B}$ does, and since event $\mathcal{B}$ occurs with strictly positive probability, we have that $\mathbb{P}(\mathcal{A}) \ge \mathbb{P}(\mathcal{B}) = c_t > 0$. Where we have labeled $\mathbb{P}(\mathcal{B})$ as $c_t$ for later convenience. □

**Appendix B. A Note on Assumption 4.2.** Here we note briefly how assumption 4.2 can be proven for arbitrary probability. Let $\mathbb{E}[\|\nabla_x^C\{F\}(\mathbf{x}) - \nabla_x^C\{G\}(\mathbf{x})\|^2] := \sigma^2$, $\mathbb{E}[\nabla_x^C\{F\}(\mathbf{x})] := \mu_F$ and $\nabla_x^C\{G\}(\mathbf{x}) := \mu_G$. For $b \ge 0$ we have $\mathbb{P}(a\|\mu_G\|^2 - \mu_F^\top \mu_G \ge b) \le \frac{\sigma^2}{b^2+\sigma^2} \le 1$ by Cantelli's concentration inequality. Therefore, $\mathbb{P}(\mu_F^\top \mu_G \ge a\|\mu_G\|^2 - b) \ge 1 - \frac{\sigma^2}{b^2+\sigma^2}$. Since $\sigma^2$ is finite, by picking $b^2 \ge \sigma^2(\frac{p}{1-p})$ we can always set $\mathbb{P}(\mu_F^\top \mu_G \ge a\|\mu_G\|^2 - b) \ge p$ for any $p$.

**Appendix C. Convergence of SMGDA in the Convex and Strongly Convex Cases.** Here we show omitted proofs of the convergence of SMGDA in the convex case and in the case where there is at least one strongly convex objective.
Let $\mathbf{x}^* = \operatorname*{argmin}_{\mathbf{s} \in P} \|\mathbf{s} - \mathbf{x}\|_2^2$, e.g. the projection of $\mathbf{x}$ to the Pareto set, then we have the following three lemmas.

LEMMA C.1. Given a set of functions $G$ with at least one $g_i m_i$-strongly convex for $i \in \{1, ..., k\}$. We then have the relation

$$682 \quad (C.1) \qquad \langle \mathbb{E}[\nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta))] - \nabla_x^C\{G\}(\mathbf{x}), \mathbf{x} - \mathbf{x}^* \rangle \ge 0$$

Proof. From the definition of convexity

$$684 \qquad \langle \mathbb{E}[\nabla_x^C\{F\}(\mathbf{x}, \mathbf{W}(\theta)) - \nabla_x^C\{G\}(\mathbf{x})], \mathbf{x} - \mathbf{x}^* \rangle. \qquad 685$$

Adding and subtracting $\sum_i \alpha_i^*(\mathbf{x}, \mathbf{W}(\theta)) \nabla g_i(\mathbf{x})$ we see that

$$687 \qquad = \langle \mathbb{E}[\sum_i \alpha_i^*(\mathbf{x}, \mathbf{W}(\theta)) \nabla f(\mathbf{x}, \mathbf{W}(\theta)) - \alpha_i^*(\mathbf{x}, \mathbf{W}(\theta)) \nabla g_i(\mathbf{x}) + (\alpha_i^*(\mathbf{x}, \mathbf{W}(\theta)) - \alpha_i^*(\mathbf{x})) \nabla g_i(\mathbf{x})], \mathbf{x} - \mathbf{x}^* \rangle,$$

$$688 \qquad = \langle \mathbb{E}[\sum_i (\alpha_i^*(\mathbf{x}, \mathbf{W}(\theta)) - \alpha_i^*(\mathbf{x})) \nabla g_i(\mathbf{x})], \mathbf{x} - \mathbf{x}^* \rangle. \qquad 689$$

Using the (strong) convexity of $g_i$ and optimality of $\alpha^*(\mathbf{x})$ we have the result,

$$691 \qquad \langle \mathbb{E}[\sum_i (\alpha_i^*(\mathbf{x}, \mathbf{W}(\theta)) - \alpha_i^*(\mathbf{x})) \nabla g_i(\mathbf{x})], \mathbf{x} - \mathbf{x}^* \rangle \ge 0. \qquad \square \\ 692$$

LEMMA C.2. Given a set of functions $G$ with at least one $g_i$ strongly convex for $i \in \{1, ..., k\}$. We have the relation

$$695 \quad (C.2) \qquad \langle \nabla_x^C\{G\}(\mathbf{x}), \mathbf{x} - \mathbf{x}^* \rangle \ge \sum_i \alpha_i^*(\mathbf{x}) \frac{m_i}{2} \|\mathbf{x} - \mathbf{x}^*\|^2$$

Proof.

$$696 \qquad \sum_i \alpha_i^*(\mathbf{x}) \langle \nabla g_i(\mathbf{x}), \mathbf{x} - \mathbf{x}^* \rangle$$
699 Using the definition of strong convexity,

$$700 \qquad \geq \sum_{i} \alpha_{i}^{*}(\mathbf{x})(g_{i}(\mathbf{x}) - g_{i}(\mathbf{x}^{*}) + \frac{m_{i}}{2} \|\mathbf{x} - \mathbf{x}^{*}\|_{2}^{2}) \qquad 701$$

702 Using the optimality of $\sum_i \alpha_i^*(\mathbf{x}^*)g_i(\mathbf{x}^*)$ we get the result,

$$703 \qquad \geq \sum_{i} \frac{\alpha_{i} m_{i}}{2} \|\mathbf{x} - \mathbf{x}^{*}\|_{2}^{2}. \qquad \square \qquad 704$$

705 LEMMA C.3. If all functions in the set G are convex, $\Phi_G(\mathbf{x})$ is also convex.

706 Proof. If $\Phi_G(\mathbf{x})$ is convex, one has the identity $\langle \nabla \Phi_G(\mathbf{x}) - \nabla \Phi_G(\mathbf{y}), \mathbf{x} - \mathbf{y} \rangle \geq 0$.

707 Starting from the definition of convexity and exploiting two reformulations of the inequality,

$$709 \qquad \langle \nabla \Phi_G(\mathbf{x}) - \nabla \Phi_G(\mathbf{y}), \mathbf{x} - \mathbf{y} \rangle = \qquad \geq \max \left( \sum_i \langle (\alpha_i^*(\mathbf{x}) - \alpha_i^*(\mathbf{y}))\nabla g_i(\mathbf{y}), \mathbf{x} - \mathbf{y} \rangle, \sum_i \langle (\alpha_i^*(\mathbf{x}) - \alpha_i^*(\mathbf{y}))\nabla g_i(\mathbf{x}), \mathbf{x} - \mathbf{y} \rangle \right). \qquad \square \qquad 710 \qquad 711$$

712 By the convexity of $g_i$ and the optimality of $\alpha^*$ we have

$$713 \qquad \langle \nabla \Phi_G(\mathbf{x}) - \nabla \Phi_G(\mathbf{y}), \mathbf{x} - \mathbf{y} \rangle \geq 0 \qquad \square \qquad 714$$

715 We can now show that if all functions in G are convex, with at least one strongly convex, that the SMGD algorithm converges to the Pareto front despite the bias.

717 THEOREM C.4. Let $g_i(\mathbf{x}) = \mathbb{E}[f_i(\mathbf{x}, \mathbf{W}(\theta))]$ be a collection of convex functions with bounded variance $\mathbb{V}(f(\mathbf{x}, \mathbf{W}(\theta))) \leq M$ and at least one $g_i$ strongly convex. Let $\mathbf{x}_t^* = \underset{s \in P-1}{\operatorname{argmin}} \|s - \mathbf{x}_t\|^2$ be the projection of $\mathbf{x}$ to the Pareto set, define $d_t^2 = \|\mathbf{x}_t - \mathbf{x}_t^*\|^2$ and a sequence $\varepsilon_t$ such that $\varepsilon_t \to 0$ as $t \to \infty$ and $\sum_t \varepsilon_t^2 \to 0$, then,

$$721 \qquad d_t^2 \xrightarrow{t} 0.$$

722 Proof. Using $\langle \nabla_x^C \{F\}(\mathbf{x}, \mathbf{W}(\theta)) - \nabla_x^C \{G\}(\mathbf{x}), d_t \rangle$ as $\langle B_t, d_t \rangle$. Starting from the
723 definition of $d_t$,

$$724 \\ 725 \qquad d_{t+1}^2 \le \|\mathbf{x}_t - \varepsilon_t \nabla_x^C \{F\}(\mathbf{x}_t) - \mathbf{x}_t^*\|_2^2.$$

726 Expanding, and adding and subtracting $\nabla_x^C\{G\}(\mathbf{x}_t)$,

$$727 \qquad \le d_t^2 - 2\varepsilon_t \langle \nabla_x^C \{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t), d_t \rangle + \varepsilon_t^2 \|\nabla_x^C \{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t)\|_2^2$$

$$728 \\ 729 \qquad \le d_t^2 - 2\varepsilon_t \langle \nabla_x^C \{G\}(\mathbf{x}_t), d_t \rangle + -2\varepsilon_t \langle B_t, d_t \rangle + \varepsilon_t^2 \|\nabla_x^C \{F\}(\mathbf{x}_t, \mathbf{W}(\theta)_t)\|_2^2$$

730 Using lemmas C.2 and C.1,

$$731 \\ 732 \qquad \mathbb{E}[d_{t+1}^2] \le \mathbb{E}[d_t^2] - 2\varepsilon_t \frac{\bar{m}}{2} \mathbb{E}[d_t^2] + \varepsilon_t^2 M.$$

733 Setting $\pi_t := \prod_{i=0}^t (1 - \varepsilon_t \bar{m})$,

$$734 \\ 735 \qquad \mathbb{E}[d_{t+1}^2] \le \pi_t d_0^2 + \sum_s \frac{\pi_t}{\pi_s} \varepsilon_s^2 M.$$

736 Requiring $\varepsilon_t \le \bar{m} \ \forall t$ and taking the limit as $t \to \infty$ gives

$$737 \\ 738 \qquad \mathbb{E}[d_{t+1}^2] \to 0. \qquad \square$$

<figure><img src="image_5.png" /></figure>

739 We can also prove convergence in the convex case (with no objective function being strongly convex).

741 THEOREM C.5. Let $g_i(\mathbf{x}) = \mathbb{E}[f_i(\mathbf{x}, \mathbf{W}(\theta))]$ be a collection of convex functions

742 with bounded variance $\mathbb{V}(f(\mathbf{x}, \mathbf{W}(\theta))) \le M$. With $\mathbf{x}_t^* = \underset{s \in P}{\operatorname{argmin}} \|s - \mathbf{x}_t\|^2$ be the

743 projection of $\mathbf{x}$ to the Pareto set, define $d_t^2 = \|\mathbf{x}_t - \mathbf{x}_t^*\|^2$ and a sequence $\varepsilon_t$ such that

744 $\varepsilon_t \to 0$ as $t \to \infty$ and $\sum_t \varepsilon_t^2 \to 0$, then,

$$745 \qquad (C.3) \qquad \mathbb{E}[\Phi_G(\mathbf{x}_t) - \Phi_G(\mathbf{x}_t^*)] \to 0 \text{ almost surely.}$$
746 Proof. Starting from the definition of $d_t$,

$$747 \qquad d_{t+1}^2 \le d_t^2 - 2\epsilon_t \langle \nabla \Phi_F(\mathbf{x}_t), \mathbf{x}_{t+1} - \mathbf{x}_t \rangle + \epsilon_t^2 \|\nabla \Phi_F(\mathbf{x}_t)\|_2^2. \qquad 748$$

749 Using $\langle \nabla_x^C \{F\}(\mathbf{x}, \mathbf{W}(\theta)) - \nabla_x^C \{G\}(\mathbf{x}), d_t \rangle$ as $\langle B_t, d_t \rangle$, taking the expectation, and
750 simplifying yields

$$751 \qquad \mathbb{E}[d_{t+1}^2] \le \mathbb{E}[d_t^2] - 2\epsilon_t \mathbb{E}\langle \nabla_x^C \{G\}(\mathbf{x}_t), d_t \rangle - 2\epsilon_t \langle B_t, d_t \rangle + \epsilon_t^2 M^2. \qquad 752$$

753 Using the convexity of $\Phi(\cdot, g)$,

$$754 \qquad \le \mathbb{E}[d_t^2] - 2\epsilon_t (\Phi_G(\mathbf{x}_t) - \Phi_G(\mathbf{x}_t^*)) + \epsilon_t^2 M^2. \qquad 755$$

756 Rearranging and summing we have,

$$757 \qquad \sum_{t=1}^{T} \epsilon_t \mathbb{E}[\Phi_G(\mathbf{x}_t) - \Phi_G(\mathbf{x}_t^*)] \le \sum_{t=1}^{T} \frac{\mathbb{E}[d_t^2 - d_{t+1}^2]}{2} + \frac{M^2}{2} \sum_{t=1}^{T} \epsilon_t^2. \qquad 758$$

759 Dividing by $E_T = \sum_{t=1}^T \epsilon_t$, we have the relation,

$$760 \qquad \sum_{t=1}^{T} \frac{\epsilon_t}{E_T} \mathbb{E}[\Phi_G(\mathbf{x}_t) - \Phi_G(\mathbf{x}_t^*)] \le \frac{d_0^2 - d_T^2}{2E_T} + \frac{M^2}{2E_T} \sum_{t=1}^{T} \epsilon_t^2. \qquad 761$$

761 Treating $\frac{\epsilon_t}{E_T}$ as a probability measure, we have by Markov's inequality

$$762 \qquad \mathbb{P}_{E_T}(\mathbb{E}[\Phi_G(\mathbf{x}_t) - \Phi_G(\mathbf{x}_t^*)] \ge \epsilon) \le \frac{d_0^2 - \mathbb{E}[d_T^2]}{2\epsilon E_T} + \frac{M^2}{2\epsilon E_T} \sum_{t=1}^{T} \epsilon_t^2 \qquad 763$$

763 The left hand side decreases to zero as $T \to \infty$, gives the result,

$$764 \qquad \mathbb{P}_{E_T}(\mathbb{E}[\Phi_G(\mathbf{x}_t) - \Phi_G(\mathbf{x}_t^*)] \ge \epsilon) \to 0. \qquad \square$$