# Fast Multiobjective Gradient Methods with Nesterov Acceleration via Inertial
Gradient-like Systems

Konstantin Sonntag  Department of Mathematics, Paderborn University, Germany
Sebastian Peitz  Department of Computer Science, Paderborn University, Germany

###### Abstract

We derive efficient algorithms to compute weakly Pareto optimal solutions for
smooth, convex and unconstrained multiobjective optimization problems in
general Hilbert spaces. To this end, we define a novel inertial gradient-like
dynamical system in the multiobjective setting, whose trajectories converge
weakly to Pareto optimal solutions. Discretization of this system yields an
inertial multiobjective algorithm which generates sequences that converge
weakly to Pareto optimal solutions. We employ Nesterov acceleration to define
an algorithm with an improved convergence rate compared to the plain
multiobjective steepest descent method (Algorithm 1). A further improvement in
terms of efficiency is achieved by avoiding the solution of a quadratic
subproblem to compute a common step direction for all objective functions,
which is usually required in first order methods. Using a different
discretization of our inertial gradient-like dynamical system, we obtain an
accelerated multiobjective gradient method that does not require the solution
of a subproblem in each step (Algorithm 2). While this algorithm does not
converge in general, it yields good results on test problems while being
faster than standard steepest descent by two to three orders of magnitude.

##  1 Introduction

In many applications in industry, economics, medicine or transport, optimizing
several criteria is of interest. In the latter, one wants to reach a
destination as fast as possible with minimal power consumption. Drug
development aims for maximizing efficacy while minimizing side effects. Even
these elementary examples share an inherent feature. The different criteria
one seeks to optimize are in general contradictory. There is no design choice
that is best for all criteria simultaneously. This insight shifts the focus
from finding a single optimal solution to a set of optimal compromises – the
Pareto set. Given the Pareto set, a decision maker can select an optimal
compromise according to their preferences. In this paper, we derive efficient
gradient-based algorithms to compute elements of the Pareto set. Formally, a
problem involving multiple criteria can be formalized as an unconstrained
multiobjective optimization problem

$\displaystyle\min_{x\in\mathcal{H}}\left[\begin{array}[]{c}f_{1}(x)\\\
\vdots\\\ f_{m}(x)\end{array}\right],$ (4)

where $f_{i}:\mathcal{H}\to\mathbb{R}$ for $i=1,\dots,m$ are the objective
functions describing the different criteria. Popular approaches to tackle this
problem in the differentiable case are first order methods which exploit the
smooth structure of the problem while not being computationally demanding
compared to higher order methods involving exact or approximated Hessians.

While in single objective optimization, _accelerated_ first order methods are
very popular, these methods are not studied sufficiently from a theoretical
point of view in the multiobjective setting. A fruitful approach to analyze
accelerated gradient methods is to interpret them as discretizations of
suitable gradient-like dynamical systems [1]. The analysis of the continuous
dynamics is often easier and can later on be transferred to the discrete
setting. This perspective is until now not fully taken advantage of in the
area of multiobjective optimization. In this paper, we utilize this approach
to derive accelerated gradient methods for multiobjective optimization. To
this end, we define and analyze the following novel dynamical gradient-like
system

$\displaystyle\ddot{x}(t)+\alpha\dot{x}(t)+\operatorname*{proj}_{C(x(t))}{(-\ddot{x}(t))}=0,$
(IMOG’)

with $\alpha>0$, $C(x)\coloneqq\operatorname{co}\left(\\{\nabla
f_{i}(x)\,:\,i=1,\dots,m\\}\right)$, where $\operatorname{co}(\cdot)$ denotes
the convex hull and $\operatorname*{proj}_{C(x(t))}{(-\ddot{x}(t))}$ is the
projection of $-\ddot{x}(t)$ onto the convex set $C(x(t))$. The system (IMOG’)
is an _inertial multiobjective gradient-like system_. We choose the
designation (IMOG’) to emphasize its relation to the system

$\displaystyle
m\ddot{x}(t)+\gamma\dot{x}(t)+\operatorname*{proj}_{C(x(t))}{(0)}=0,$ (IMOG)

with $m,\gamma>0$, which was discussed in [2]. In the single objective setting
($m=1$), both (IMOG’) and (IMOG) reduce to the _heavy ball with friction
system_

$\displaystyle m\ddot{x}(t)+\gamma\dot{x}(t)+\nabla f(x(t))=0,$ (HBF)

which is well studied for different types of objective functions $f$, see e.g.
[3, 4, 5].

We discretize (IMOG’) to obtain an iterative scheme of the form

$\displaystyle
x^{k+1}=x^{k}+a(x^{k}-x^{k-1})-s\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k}),$

with appropriately chosen coefficients $a,s>0$ and
$\theta^{k}\in\mathbb{R}^{m}$. This scheme can be interpreted as an inertial
gradient method for (4). We show that it shares many properties with its
continuous counterpart and that iterates defined by this algorithm converge
weakly to Pareto critical points. To the best of our knowledge this is the
first multiobjective method involving a constant momentum term with guaranteed
convergence to Pareto critical solutions.

In a further step, we introduce time-dependent friction and informally define
the following _multiobjective gradient-like system with asymptotically
vanishing damping_

$\displaystyle\ddot{x}(t)+\frac{\alpha}{t}\dot{x}(t)+\operatorname*{proj}_{C(x(t))}{(-\ddot{x}(t))}=0.$
(MAVD)

In the single objective setting this system simplifies to the following
inertial system with _asymptotically vanishing damping_

$\displaystyle\ddot{x}(t)+\frac{\alpha}{t}\dot{x}(t)+\nabla f(x(t))=0.$ (AVD)

It is well-known that (AVD) is naturally linked with Nesterov’s accelerted
gradient method [1, 6, 7, 8]. Discretizing the dynamical system (MAVD), and
using our knowledge about (IMOG’), we derive an accelerated gradient method
for multiobjective optimization that takes the form

$\displaystyle
x^{k+1}=x^{k}+\frac{k-1}{k+2}(x^{k}-x^{k-1})-s\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k}),$

with appropriately chosen coefficients $s>0$ and
$\theta^{k}\in\mathbb{R}^{m}$. In a recent preprint, Tanabe, Fukuda and
Yamashita derive an accelerated proximal gradient method for multiobjective
optimization using the concept of merit functions [9]. We show that the method
we derive from the differential equation (MAVD) achieves the same convergence
rate of order $\mathcal{O}(k^{-2})$ for the function values, measured with a
so-called merit function.

This remainder of the paper is organized according to the novelties described
above. After introducing some basic definitions and notation in Section 2, we
prove that solutions to the system (IMOG’) exist in finite-dimensional Hilbert
spaces in Section 3, and that they converge to Pareto critical points in
Section 4. Based on that, we derive a discrete optimization algorithm from an
explicit discretization on (IMOG’) and show that the iterates defined by this
method converge weakly to Pareto critical points, in Section 5. We then
introduce Nesterov Acceleration and prove an improved convergence result in
Section 6, and a further improvement in numerical efficiency is discussed in
Section 7. These two central algorithms are summarized in Algorithms 1 and 2
in the respective sections. We compare these on convex and non-convex example
problems in Section 8, and conclude our findings and list future research
projects in Section 9.

##  2 Background

###  2.1 Notations

Throughout this paper, $\mathcal{H}$ is a real Hilbert space with inner
product $\langle\cdot,\cdot\rangle$ and induced norm $\lVert\cdot\rVert$. The
set $\Delta_{m}\coloneqq\left\\{\alpha\in\mathbb{R}^{m}\,\,:\,\,\alpha\geq
0,\,\,\text{and}\,\,\sum_{i=1}^{m}\alpha_{i}=1\right\\}$ is the unit simplex.
We denote the open ball with radius $\delta>0$ and center $x$ by
$B_{\delta}(x)\coloneqq\left\\{y\in\mathcal{H}\,\,:\,\,\lVert
y-x\rVert<\delta\right\\}$. The closed ball with radius $\delta>0$ and center
$x$ is denoted by $\overline{B_{\delta}(x)}$. For a set of vectors
$\left\\{\xi_{1},\cdots,\xi_{k}\right\\}\subseteq\mathcal{H}$ we denote the
convex hull of these vectors by
$\operatorname{co}(\left\\{\xi_{1},\cdots,\xi_{k}\right\\})\coloneqq\left\\{\sum_{i=1}^{m}\alpha_{i}\xi_{i}\,\,:\,\,\alpha\in\Delta_{m}\right\\}$.
For a closed convex set $C\subseteq\mathcal{H}$ the projection of a vector
$x\in\mathcal{H}$ onto $C$ is
$\operatorname*{proj}_{C}(x)\coloneqq\operatorname*{arg\,min}_{y\in\mathcal{H}}\lVert
y-x\rVert^{2}$. For two vectors $x,y\in\mathbb{R}^{k}$, we define the partial
order $x\leq y:\Leftrightarrow x_{i}\leq y_{i}$ for all $i=1,\dots,k$. We
define $\geq$, $<$, $>$ on $\mathbb{R}^{k}$ analogously. When we treat
dynamical systems, $t\in\mathbb{R}$ and $x\in\mathcal{H}$ are the time and
state variable, respectively. We denote trajectories in $\mathcal{H}$ with
$t\mapsto x(t)$ with first derivative $\dot{x}(t)$ and second derivative
$\ddot{x}(t)$.

###  2.2 Multiobjective Optimization

Consider the unconstrained multiobjective optimization problem

$\displaystyle\min_{x\in\mathcal{H}}\left[\begin{array}[]{c}f_{1}(x)\\\
\vdots\\\ f_{m}(x)\end{array}\right],$ (8)

with at least once differentiable objective functions
$f_{i}:\mathcal{H}\to\mathbb{R}$ for $i=1,\dots,m$. The definitions in this
subsection are aligned with [10].

######  Definition 2.1 ([10], pp.10-20).

Consider the optimization problem (8).

  1. i)

A point $x^{*}\in\mathcal{H}$ is Pareto optimal if there does not exist
another point $x\in\mathcal{H}$ such that $f_{i}(x)\leq f_{i}(x^{*})$ for all
$i=1,\dots,m,$ and $f_{j}(x)<f_{j}(x^{*})$ for at least one index $j$. The set
of all Pareto optimal points is the Pareto set, which we denote with $P$.

  2. ii)

A point $x^{*}\in\mathcal{H}$ is locally Pareto optimal if there exists
$\delta>0$ such that $x^{*}$ is Pareto optimal in $B_{\delta}(x^{*})$.

  3. iii)

A point $x^{*}\in\mathcal{H}$ is weakly Pareto optimal if there does not exist
another vector $x\in\mathcal{H}$ such that $f_{i}(x)<f_{i}(x^{*})$ for all
$i=1,\dots,m$.

  4. iv)

A point $x^{*}\in\mathcal{H}$ is locally weakly Pareto optimal if there exists
$\delta>0$ such that $x^{*}$ is weakly Pareto optimal in $B_{\delta}(x^{*})$.

In this paper, we treat convex MOPs. Therefore, the objective functions
$f_{i}$ are convex for all $i=1,\dots,m$. In this setting, every locally
(weakly) Pareto optimal point is also (weakly) Pareto optimal. For
unconstrained MOPs, the so-called Karush-Kuhn-Tucker conditions can be written
as follows.

######  Definition 2.2.

A point $x^{*}\in\mathcal{H}$ satisfies the Karush-Kuhn-Tucker conditions if
there exists $\alpha\in\Delta_{m}$ such that $\sum_{i=1}^{m}\alpha_{i}\nabla
f_{i}(x^{*})=0$. If $x^{*}$ satisfies the Karush-Kuhn-Tucker conditions, we
call it Pareto critical.

The condition $0\in\operatorname{co}\left(\left\\{\nabla
f_{i}(x^{*})\,\,:\,\,i=1,\dots,m\right\\}\right)$ is equivalent to the Karush-
Kuhn-Tucker conditions. Analogously to the single objective setting,
criticality of a point is a necessary condition for optimality. In the convex
setting, the KKT conditions are also sufficient conditions for weak Pareto
optimality. If we denote the Pareto set by $P$, the weak Pareto set by $P_{w}$
and the Pareto critical set by $P_{c}$ in the setting of smooth and convex
multiobjective optimization, we observe the relation

$\displaystyle P\subset P_{w}=P_{c}.$

###  2.3 Accelerated Methods for Multiobjective Optimization

Accelerated methods for multiobjective optimization are not sufficiently
discussed from a theoretical point of view in the literature yet. In [11] El
Moudden and El Moutasim propose an accelerated method for multiobjective
optimization which incorporates the multiobjective descent direction by Fliege
[12] and the same acceleration scheme as in Nesterov’s accelerated method
[13]. El Moudden and El Moutasim prove a convergence rate of the function
values with rate $\mathcal{O}(k^{-2})$. Their proof relies on the restrictive
assumption that the Lagrange multipliers of the quadratic subproblem, that is
used to compute the step direction in every iteration, remain fixed from a
certain point on. Under this assumption, the method simplifies to Nesterov’s
method for single objective optimization problems applied to a weighted sum of
the objective functions with fixed weights. Only recently, Tanabe, Fukuda and
Yamashita derived an accelerated proximal gradient method for multiobjective
optimization problems in [9]. They developed their method using the concept of
merit functions (see Subsection 2.5) and show that the function values
converge with rate $\mathcal{O}(k^{-2})$ without additional assumptions on the
Lagrange multipliers.

###  2.4 Dynamical Systems linked to Multiobjective Optimization

In [14] Smale presents the idea of treating multiobjective optimization
problem with a continuous time perspective that is motivated from an
economical point of view using utility functions in a multi-agent framework.
The simplest dynamical system for multiobjective optimization problems is the
_multiobjective gradient system_

$\displaystyle\dot{x}(t)+\operatorname*{proj}_{C(x(t))}(0)=0,$ (MOG)

where $C(x(t))=\operatorname{co}\left(\\{\nabla
f_{i}(x(t))\,:\,i=1,\dots,m\\}\right)$. This system is already treated in [15]
and in addition by Cornet in [16, 17, 18] as a dynamical system for resource
allocation from an economical point of view. In [19, 20] the system (MOG) gets
introduced as a tool for multiobjective optimization. The system (MOG) can
also be seen as a continuous version of the multiobjective steepest descent
method by Fliege [12]. In the single objective setting ($m=1$), the system
(MOG) simplifies to the _steepest descent dynamical system_ $\dot{x}(t)+\nabla
f(x(t))=0$. Generalizations of (MOG) are treated in [20, 21]. In [20] Attouch
and Goudou discuss a dynamical system for constrained minimization and in [21]
Attouch, Garrigos and Goudou present a differential inclusion for constrained
nonsmooth optimization.
In [2], Attouch and Garrigos introduce inertia in the system (MOG) and define
the following _inertial multiobjective gradient-like dynamical system_

$\displaystyle
m\ddot{x}(t)+\gamma\dot{x}(t)+\operatorname*{proj}_{C(x(t))}(0)=0.$ (IMOG)

Trajectories of (IMOG) converge weakly to Pareto optimal solutions given
$\gamma^{2}>mL$, where $L$ is a common Lipschitz constant for the gradients of
the objective functions.

###  2.5 Merit Functions

A merit function associated with an optimization problem is a function that
returns zero at an optimal solution and which is strictly positive otherwise.
An overview on merit functions used in multiobjective optimization is given in
[22]. In our proofs we use the merit function

$\displaystyle
u_{0}(x)\coloneqq\sup_{z\in\mathcal{H}}\min_{i=1,\dots,m}f_{i}(x)-f_{i}(z),$
(9)

which satisfies the following statement.

######  Theorem 2.3 ([22] Theorem 3.1).

It holds that $u_{0}(x)\geq 0$ for all $x\in\mathcal{H}$. Moreover,
$x\in\mathcal{H}$ is weakly Pareto optimal for (8), if and only if
$u_{0}(x)=0$.

Additionally, $u_{0}(x)$ is lower semicontinuous. Therefore, if
$(x^{k})_{k\geq 0}$ is a sequence with $u_{0}(x^{k})\to 0$, every cluster
point of $(x^{k})_{k\geq 0}$ is weakly Pareto optimal. This motivates the
usage of $u_{0}(x)$ as a measure of complexity for multiobjective optimization
methods. The function $u_{0}(x)$ is not the only merit function for
multiobjective optimization problems, see also [23, 24, 25] and further
references in [22].

##  3 Global Existence in Finite Dimensions

In this section, we show that solutions exist for the Cauchy problem related
to (IMOG’), i.e

$\displaystyle\left\lvert\begin{array}[]{l}\ddot{x}(t)+\alpha\dot{x}(t)+\operatorname*{proj}_{C(x(t))}{(-\ddot{x}(t))}=0,\\\
\\\ x(0)=x_{0},\quad\dot{x}(0)=v_{0},\end{array}\right.$ (13)

with initial data $x_{0},v_{0}\in\mathcal{H}$. To this end, we show that for
this system solutions exist if there exists a solution for a first order
differential inclusion

$\displaystyle(\dot{u},\dot{v})\in F(u,v),$

with a set-valued mapping
$F:\mathcal{H}\times\mathcal{H}\rightrightarrows\mathcal{H}\times\mathcal{H}$.
Then, we use an existence theorem for differential inclusions from [26]. Our
argument works only in finite-dimensional Hilbert spaces. Thus, we assume
$\dim(\mathcal{H})<+\infty$ from here on. In our context, the following set-
valued mapping is of interest:

$\displaystyle
F:\mathcal{H}\times\mathcal{H}\rightrightarrows\mathcal{H}\times\mathcal{H},\quad(u,v)\mapsto\\{v\\}\times\left(\big{(}-\operatorname*{arg\,min}_{z\in
C(u)}\langle z,v\rangle\big{)}-\alpha v\right)$ (14)

As stated above, $C(u)\coloneqq\operatorname{co}\left(\\{\nabla
f_{i}(u)\,:\,i=1,\dots,m\\}\right)$. We can show that (13) has a solution if
the differential inclusion

$\displaystyle\left\lvert\begin{array}[]{l}(\dot{u}(t),\dot{v}(t))\in
F(u,v),\\\ \\\ (u(0),v(0))=(u_{0},v_{0}),\end{array}\right.$ (18)

with appropriate initial data $u_{0},v_{0}$ has a solution.

###  3.1 Existence of Solutions to (DI)

To show that there exist solutions to (18), we investigate the set-valued
mapping $(u,v)\rightrightarrows F(u,v)$ defined in (14). The basic definitions
for set-valued mappings used in this subsection can be found in [26].

######  Proposition 3.1.

For all $(u,v)\in\mathcal{H}\times\mathcal{H}$,
$F(u,v)\subset\mathcal{H}\times\mathcal{H}$ is convex and compact.

###### Proof.

The statement follows directly from the definition. ∎

To use an existence theorem from [26], we need to show that
$(u,v)\rightrightarrows F(u,v)$ is upper semicontinuous. Showing this is
elementary. We omit the full proof here but sketch a possible way to prove
this result.

######  Lemma 3.2.

Let $C(u)\coloneqq\operatorname{co}(\\{c_{i}(u)\,:\,i=1,\dots,m\\})$ with
$c_{i}:\mathcal{H}\to\mathcal{H},\,u\to c_{i}(u)$ continuous for
$i=1,\dots,m$. Let $(u_{0},v_{0})\in\mathcal{H}\times\mathcal{H}$ be fixed.
Then, for all $\varepsilon>0$ there exists an $\delta>0$ such that for all
$(u,v)\in\mathcal{H}\times\mathcal{H}$ with $\lVert u-u_{0}\rVert<\delta$ and
$\lVert v-v_{0}\rVert<\delta$ and for all $z\in\operatorname*{arg\,min}_{z\in
C(u)}\langle z,v\rangle$ there exists
$z_{0}\in\operatorname*{arg\,min}_{z_{0}\in C(u_{0})}\langle
z_{0},v_{0}\rangle$ with $\lVert z-z_{0}\rVert<\varepsilon$.

###### Proof.

The proof follows straightforward. ∎

######  Proposition 3.3.

The set-valued mapping $(u,v)\rightrightarrows F(u,v)$ is upper
semicontinuous.

###### Proof.

Using Lemma 3.2 we can show in a straightforward manner

$\displaystyle F((u_{0},v_{0})+B_{\delta}(0))\subset
F(u_{0},v_{0})+B_{\varepsilon}(0)$

using only continuity arguments. Then, the statement follows by the fact that
$(u,v)\rightrightarrows F(u,v)$ is locally compact. ∎

######  Proposition 3.4.

Let $\mathcal{H}$ have finite dimension. Then, the mapping

$\displaystyle\phi:\mathcal{H}\times\mathcal{H}\to\mathcal{H}\times\mathcal{H},\quad(u,v)\mapsto\left(v,\,\operatorname*{proj}_{F(u,v)}(0)\right),$

is locally compact.

###### Proof.

If $\dim(\mathcal{H})<+\infty$ the proof follows easily since all images
$F(u,v)$ are compact and depend on $(u,v)$ in a well behaved manner. On the
other hand, from $\phi$ being locally compact, we get that $v\mapsto v$ is
locally compact which is equivalent to $\mathcal{H}$ being finite-dimensional.
∎

The following existence theorem from [26] is applicable in our setting.

######  Theorem 3.5 ([26], p. 98, Theorem 3).

Let $X$ be a Hilbert space, $\Omega\subset\mathbb{R}\times X$ be an open
subset containing $(0,x_{0})$. Let $F$ be an upper semicontinuous map from
$\Omega$ into the non-empty closed convex subsets of $X$. We assume that
$(t,x)\mapsto\operatorname*{proj}_{F(t,x)}(0)$ is locally compact. Then, there
exists $T>0$ and an absolutely continuous function $x(\cdot)$ defined on
$[0,T]$ which is a solution to the differential inclusion

$\displaystyle\dot{x}(t)\in F(t,x(t)),\quad x(0)=x_{0}.$

We are finally in the position to state an existence theorem for (18).

######  Theorem 3.6.

Assume $\mathcal{H}$ is finite-dimensional and that the gradients of the
objective function $\nabla f_{i}$ are globally Lipschitz continuous. Then, for
all $(u_{0},v_{0})\in\mathcal{H}\times\mathcal{H}$ there exists $T>0$ and an
absolutely continuous function $(u(\cdot),v(\cdot))$ defined on $[0,T]$ which
is a solution to the differential inclusion (18).

###### Proof.

The proof follows immediately from Propositions 3.3 and 3.4 which show that
the set-valued mapping satisfies all conditions required for Theorem 3.5. ∎

In the following, we show that under additional conditions on the objective
functions $f_{i}$ there exist solutions defined on $[0,+\infty)$. The
extension of the solution works using a standard argument. We show that the
solutions to (18) remain bounded. Then, we use Zorn’s Lemma to retrieve a
contradiction if there is a maximal solution that is not defined on
$[0,+\infty)$.

######  Theorem 3.7.

Assume $\mathcal{H}$ is finite-dimensional and that the gradients of the
objective function $\nabla f_{i}$ are globally Lipschitz continuous. Then, for
all $(u_{0},v_{0})\in\mathcal{H}\times\mathcal{H}$ there exists an absolutely
continuous function $(u(\cdot),v(\cdot))$ defined on $[0,+\infty)$ which is a
solution to the differential inclusion (18).

###### Proof.

Theorem 3.6 guarantees the existence of solutions defined on $[0,T)$ for some
$T\geq 0$. Using the domain of definition, we can define a partial order on
the set of solutions to the problem (18). Assuming there is no solution
defined on $[0,+\infty)$, Zorn’s Lemma guarantees the existence of a solution
$(u(\cdot),v(\cdot)):[0,T)\to\mathcal{H}\times\mathcal{H}$ with $T<+\infty$
which can not be extended. We will show that $(u(t),v(t))$ does not blow up in
finite time and therefore can be extended which contradicts the claimed
maximality.
Define

$\displaystyle
h(t)\coloneqq\lVert(u(t),v(t))-(u(0),v(0))\rVert_{\mathcal{H}\times\mathcal{H}},$

where $\lVert(x,y)\rVert_{\mathcal{H}\times\mathcal{H}}=\sqrt{\lVert
x\rVert^{2}+\lVert y\rVert^{2}}$. We show that $h(t)$ can be bounded by a
real-valued function. Using the Cauchy-Schwarz inequality, we directly get

$\displaystyle\frac{d}{dt}\frac{1}{2}h^{2}(t)$
$\displaystyle=\langle(\dot{u}(t),\dot{v}(t)),(u(t),v(t))-(u(0),v(0))\rangle$
(19) $\displaystyle\leq\lVert F(u(t),v(t))\rVert h(t).$ (20)

We next derive a bound on $\lVert F(u,v)\rVert$. The basic inequalities
between the $\ell_{1}$ and $\ell_{2}$ norm applied to $\left(\lVert
x\rVert,\lVert y\rVert\right)\in\mathbb{R}^{2}$ yield

$\displaystyle\lVert(x,y)\rVert_{\mathcal{H}\times\mathcal{H}}\leq\lVert
x\rVert+\lVert
y\rVert\leq\sqrt{2}\lVert(x,y)\rVert_{\mathcal{H}\times\mathcal{H}}.$

Therefore

$\displaystyle\lVert F(u,v)\rVert_{\mathcal{H}\times\mathcal{H}}\leq\lVert
v\rVert+\lVert\big{(}-\operatorname*{arg\,min}_{z\in C(u)}\langle
z,v\rangle\big{)}-\alpha v\rVert$ $\displaystyle\leq$
$\displaystyle(1+\alpha)\lVert
v\rVert+\max_{\theta\in\Delta_{m}}\lVert\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(u)\rVert$ $\displaystyle\leq$ $\displaystyle(1+\alpha)\lVert
v\rVert+\max_{\theta\in\Delta_{m}}\lVert\sum_{i=1}^{m}\theta_{i}\left(\nabla
f_{i}(u)-\nabla
f_{i}(0)\right)\rVert+\max_{\theta\in\Delta_{m}}\lVert\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(0)\rVert$ $\displaystyle\leq$ $\displaystyle(1+\alpha)\lVert
v\rVert+L\lVert u\rVert+\max_{i=1,\dots,m}\lVert\nabla f_{i}(0)\rVert$
$\displaystyle\leq$ $\displaystyle
c(1+\lVert(u,v)\rVert_{\mathcal{H}\times\mathcal{H}}),$

where we have chosen
$c=\sqrt{2}\max\left\\{1+\alpha,L,\max_{i=1,\dots,m}\lVert\nabla
f_{i}(0)\rVert\right\\}$. Combining the last inequality with (19), we get

$\displaystyle\frac{d}{dt}\frac{1}{2}h^{2}(t)\leq c(1+h(t))h(t),\,\,\text{for
all}\,\,t\in[0,T).$

Using a Gronwall-type argument (see Lemma A.4 and Lemma A.5 in [27]) just as
in Theorem 3.5 in [2], we know that there exists $c\in\mathbb{R}$ such that
for an arbitrary $\varepsilon>0$

$\displaystyle h(t)\leq cT\exp(cT),\,\,\text{ for all
}\,\,t\in[0,T-\varepsilon].$

Since this upper bound is independent of $t$ and $\varepsilon$, it follows
that $h\in L^{\infty}([0,T],\mathbb{R})$. Therefore, solutions to (18) do not
blow up in finite time and can be extended. This is a contradiction to the
maximality of the solution $(u(t),v(t))$. ∎

###  3.2 Existence of Solutions to (CP)

Using the findings of the previous subsection, we can proceed with the
discussion of the Cauchy problem (13). In this subsection, we show that
solutions to the differential inclusion (18) immediately give solutions to the
Cauchy problem (13).

######  Theorem 3.8.

Let $x_{0},v_{0}\in\mathcal{H}$. Assume that $(u(t),v(t))$ for
$t\in[0,+\infty)$ is a solution to (18) with $(u(0),v(0))=(x_{0},v_{0})$.
Then, it follows that $x(t)\coloneqq u(t)$ satisfies the differential equation

$\displaystyle\ddot{x}(t)+\alpha\dot{x}(t)+\operatorname*{proj}_{C(x(t))}(-\ddot{x}(t))=0,\,\,\text{for
all}\,\,t\in(0,+\infty),$

and $x(0)=x_{0}$, $\dot{x}(0)=v_{0}$.

###### Proof.

Since $(u(t),v(t))$ is a solution to (18), it holds for all $t\in(0,+\infty)$
that

$\displaystyle\dot{u}(t)$ $\displaystyle=v(t),$ $\displaystyle\dot{v}(t)$
$\displaystyle\in\operatorname*{arg\,min}_{z\in C(u(t))}\langle
z,v(t)\rangle-\alpha v(t).$

Using Lemma A.1, the second line gives $-\alpha
v(t)=\operatorname*{proj}_{C(u(t))+\dot{v}(t)}(0)$, which is equivalent to

$\displaystyle\dot{v}(t)+\alpha
v(t)+\operatorname*{proj}_{C(u(t))}(-\dot{v}(t))=0.$

Rewriting this system using $x(t)=u(t)$, $\dot{x}(t)=\dot{u}(t)=v(t)$ and
$\ddot{x}=\dot{v}(t)$ and verifying the initial conditions $x(0)=u(0)=x_{0}$
and $\dot{x}(0)=v(0)=v_{0}$ yields the desired result. ∎

Finally, we can state the full existence theorem for the Cauchy problem (13).

######  Theorem 3.9.

Assume $\mathcal{H}$ is finite-dimensional and that the gradients of the
objective function $\nabla f_{i}$ are globally Lipschitz continuous. Then, for
all $x_{0},v_{0}\in\mathcal{H}$, there exists a twice continuously
differentiable function $x(t)$ defined on $[0,+\infty)$ which is absolutely
continuous with absolutely continuous first derivative $\dot{x}(t)$, and which
is a solution to the Cauchy problem (13) with initial values $(x_{0},v_{0})$.

###### Proof.

The proof follows immediately combining Theorem 3.7 and Theorem 3.8. ∎

######  Remark 3.10.

Throughout this section, we have assumed that the gradients $\nabla f_{i}$ of
the objective functions are globally Lipschitz continuous. One can relax this
condition and only require the gradients to be Lipschitz continuous on bounded
sets if we can guarantee that the solutions remain bounded. This holds for
example if one of the objective functions $f_{i}$ is coercive.

##  4 Asymptotic Analysis of Trajectories of (IMOG’)

In this section, we omit the assumption $\dim(\mathcal{H})<+\infty$. We show
that trajectories of the differential equation (IMOG’) converge weakly to
Pareto critical points of the optimization problem (8). This follows from a
dissipative property of the system and an argument that relies on Opial’s
Lemma. We first define a so-called energy function for the system (IMOG’) that
has Lyapunov-type properties.

######  Proposition 4.1.

Let $x:[0,+\infty)\to\mathcal{H}$ be a solution to (13). For $i=1,\dots,m$
define the global energy

$\displaystyle\mathcal{E}_{i}:[0,T)\to\mathbb{R},\quad t\mapsto
f_{i}(x(t))+\frac{1}{2}\lVert\dot{x}(t)\rVert^{2}.$

Then, for all $t\in(0,+\infty)$ it holds that
$\frac{d}{dt}\mathcal{E}_{i}(t)\leq-\alpha\lVert\dot{x}(t)\rVert^{2}$.

###### Proof.

From the definition of the differential equation (IMOG’), it follows that
$-\alpha\dot{x}(t)=\operatorname*{proj}_{\operatorname{co}\left(\nabla
f_{i}(x(t))+\ddot{x}(t)\right)}(0)$. By the variational characterization of
the convex projection, we get for all $i=1,\dots,m$

$\displaystyle\langle\alpha\dot{x}(t)+\nabla
f_{i}(x(t))+\ddot{x}(t),\alpha\dot{x}(t)\rangle\leq 0,$

which immediately gives

$\displaystyle\langle\nabla
f_{i}(x(t)),\dot{x}(t)\rangle+\langle\dot{x}(t),\ddot{x}(t)\rangle\leq-\alpha\lVert\dot{x}(t)\rVert^{2}.$

Applying the chain rule to $\frac{d}{dt}\mathcal{E}_{i}(t)$ yields the desired
result. ∎

######  Proposition 4.2.

Let $x:[0,+\infty)\to H$ be a bounded solution of (13) and let further $\nabla
f_{i}$ be Lipschitz continuous on bounded sets. Then, for all $i=1,\dots,m$

  1. i)

$\lim_{t\to+\infty}\mathcal{E}_{i}(t)=\mathcal{E}_{i}^{\infty}>-\infty$

  2. ii)

$\dot{x}\in L^{2}([0,+\infty))\cap L^{\infty}([0,+\infty))$

  3. iii)

$\ddot{x}\in L^{\infty}([0,+\infty))$,
$\lim_{t\to+\infty}\lVert\dot{x}(t)\rVert=0$ and
$\lim_{t\to+\infty}f_{i}(x(t))=\mathcal{E}_{i}^{\infty}$

  4. iv)

There exists a monotonically increasing unbounded sequence $(t_{k})_{k\geq 0}$
with
$\operatorname*{proj}_{C(x(t_{k}))}(0)\to 0$ for $k\to+\infty$.

###### Proof.

i) From Proposition 4.1, we immediately get that $\mathcal{E}_{i}$ is
monotonically decreasing and therefore
$\mathcal{E}_{i}(t)\to\mathcal{E}_{i}^{\infty}$ as $t\to+\infty$. We have to
show that in fact $\mathcal{E}_{i}^{\infty}>-\infty$. Since $\nabla f_{i}$ is
bounded on bounded sets, we can conclude by the mean value theorem that
$f_{i}$ is bounded on bounded sets. Since $x(t)$ remains bounded by
assumption, we conclude that $f_{i}(x(t))$ is bounded from below, and hence

$\displaystyle\mathcal{E}_{i}^{\infty}\geq\inf_{t\geq 0}f_{i}(x(t))>-\infty.$

ii) We know that $f_{i}(x(t))$ is bounded. Then, by the definition of
$\mathcal{E}_{i}$ and the fact that $\mathcal{E}_{i}$ is monotonically
decreasing, we immediately get that $\dot{x}$ is bounded for all $t\geq 0$.
Since it is of class $C^{1}$, it follows that $\dot{x}\in
L^{\infty}([0,+\infty))$. Using Proposition 4.1 it follows that

$\displaystyle\alpha\int_{0}^{+\infty}\lVert\dot{x}(t)\rVert^{2}\,dt\leq-\int_{0}^{+\infty}\frac{d}{dt}\mathcal{E}_{i}(s)\,ds=\mathcal{E}_{i}(0)-\mathcal{E}_{i}^{\infty}<+\infty,$

and therefore $\dot{x}\in L^{2}([0,+\infty))$.

iii) Since $\dot{x}(t)$ and $\nabla f_{i}(x(t))$ remain bounded for all $t\geq
0$ it follows that
$\ddot{x}(t)=-\alpha\dot{x}(t)-\operatorname*{proj}_{C(x(t))}(-\ddot{x}(t))$
remains bounded for all $t\geq 0$. By the fact that $x$ is twice continuously
differentiable, it follows that $\ddot{x}$ is continuous and hence
$\ddot{x}\in L^{\infty}([0,+\infty))$. Then, from $\dot{x}\in
C^{1}([0,+\infty))\cap L^{2}([0,+\infty))$ and $\ddot{x}\in C([0,+\infty))\cap
L^{\infty}([0,+\infty))$ it follows that
$\lim_{t\to+\infty}\lVert\dot{x}(t)\rVert=0$. From
$\lim_{t\to+\infty}\lVert\dot{x}(t)\rVert=0$ and part i) we can immediately
conclude $\lim_{t\to+\infty}f_{i}(x(t))=\mathcal{E}_{i}^{\infty}$.

iv) Assume that the negation of statement iv) holds, namely

$\displaystyle\exists M>0\,\,\exists T\geq 0\,\,\forall t\geq
T\,\,:\,\,\lVert\operatorname*{proj}_{C(x(t))}(0)\rVert\geq M.$ (21)

Fix an arbitrary $\delta>0$. Since $\dot{x}(t)\to 0$ and $\nabla f_{i}$ is
Lipschitz continuous on a set containing $x(t)$ it follows that there exists
$T_{\delta}>T$ such that for all $t>T_{\delta}$ and $s\in[t,t+\delta)$

$\displaystyle\begin{split}\lVert\nabla f_{i}(x(s))-\nabla
f_{i}(x(t))\rVert<\frac{M}{8}\,\,\text{and}\,\,\lVert\alpha\dot{x}(s)\rVert<\frac{M}{8}.\end{split}$
(22)

Define
$v\coloneqq\operatorname*{proj}_{C(x(t))}/\lVert\operatorname*{proj}_{C(x(t))}\rVert$.
From (21) it follows that

$\displaystyle\forall\xi\in C(x(t))\,\,:\,\,\langle\xi,v\rangle\geq M.$

Combining the last statement with (22) and using the Cauchy-Schwarz inequality
we get

$\displaystyle\forall s\in[t,t+\delta)\,\,\forall\xi\in
C(x(s))\,\,:\,\,\langle\xi+\alpha\dot{x}(s),v\rangle\geq\frac{M}{4}.$

And hence

$\displaystyle\forall
s\in[t,t+\delta)\,\,:\,\,\langle-\ddot{x}(s),v\rangle\geq\frac{M}{4}.$

Using the Cauchy-Schwarz inequality again, we get

$\displaystyle\lVert\dot{x}(t)-\dot{x}(t+\delta)\rVert\geq\langle\dot{x}(t)-\dot{x}(t+\delta),v\rangle=\langle\int_{t}^{t+\delta}-\ddot{x}(s)\,ds,v\rangle$
$\displaystyle=$
$\displaystyle\int_{t}^{t+\delta}\langle-\ddot{x}(s),v\rangle\,ds\geq\int_{t}^{t+\delta}\frac{M}{4}\,ds=\frac{M\delta}{4}.$

Since we can choose $\delta$ arbitrarily large and independently from $M$,
this contradicts $\dot{x}(t)\to 0$. ∎

We will use part iv) of Proposition 4.2 to show that a weak limit point of the
trajectory $x(t)$ is Pareto critical. To this end we introduce the following
lemma that states a demiclosedness property of the set-valued mapping
$C:\mathcal{H}\rightrightarrows\mathcal{H},\,x\mapsto C(x)$.

######  Lemma 4.3 ([6] Lemma 2.4, [2] Lemma 4.10).

Assume that the objective functions $f_{i}$ are continuously differentiable.
Let $(x^{k})_{k\geq 0}$ be a sequence in $\mathcal{H}$ that converges weakly
to $x^{\infty}$, and assume there exists a sequence $(g^{k})_{k\geq 0}$ with
$g^{k}\in C(x^{k})$ that converges strongly to zero. Then, $0\in
C(x^{\infty})$ and hence $x^{\infty}$ is Pareto critical.

If we can show that the trajectories of (IMOG’) converge weakly, Proposition
4.2 together with Lemma 4.3 guarantees that the limit points are Pareto
critical. To show that the trajectories are in fact converging, we require
Opial’s Lemma.

######  Lemma 4.4 (Opial’s Lemma [28]).

Let $S\subset\mathcal{H}$ be a nonempty subset of $\mathcal{H}$ and
$x:[0,+\infty)\to\mathcal{H}$. Assume that $x(t)$ satisfies the following
conditions.

  1. i)

Every weak sequential cluster point of $x(t)$ belongs to $S$.

  2. ii)

For every $z\in S$, $\lim_{t\to+\infty}\lVert x(t)-z\rVert$ exists.

Then, $x(t)$ converges weakly to an element $x^{\infty}\in S$.

To use Opial’s Lemma, we need a suitable nonempty set $S\subset\mathcal{H}$
that we define in the following proposition.

######  Proposition 4.5.

Let $x(t)$ be a bounded solution to (13). Then, the set

$\displaystyle
S\coloneqq\left\\{z\in\mathcal{H}\,:\,f_{i}(z)\leq\mathcal{E}_{i}^{\infty}\,\,\text{for
all }i=1,\dots,m,\right\\},$ (23)

is nonempty.

###### Proof.

Part iii) of Proposition 4.2 states that
$\lim_{t\to+\infty}f_{i}(x(t))=\mathcal{E}_{i}^{\infty}$ for all
$i=1,\dots,m$. Since $x(t)$ is bounded, it possesses at least one weak
sequential cluster point $x^{\infty}$. The objective functions $f_{i}$ are
convex and continuous and therefore weakly lower semicontinuous. From this we
conclude $x^{\infty}\in S$. ∎

For the set $S$ defined in (23) and a bounded solution $x(t)$ of (13), the
first part of Opial’s Lemma is easy to obtain. It follows analogously to the
proof of Proposition 4.5 where it is shown that $S$ is nonempty. To show the
second part of Opial’s Lemma, we verify that
$h_{z}(t)\coloneqq\frac{1}{2}\lVert x(t)-z\rVert^{2}$ satisfies a differential
inequality. Then, the convergence can be deduced from the following lemma.

######  Lemma 4.6 ([5] Lemma 4.2).

Let $h\in C^{1}([0,+\infty),\mathbb{R})$ be a positive function satisfying
$\alpha\dot{h}(t)+\ddot{h}(t)\leq g(t)$ for all $t\geq 0$, with $g\in
L^{1}([0,+\infty),\mathbb{R})$ and $\alpha>0$. Then, $\lim_{t\to+\infty}h(t)$
exists.

With these ingredients, we can formulate the main convergence theorem of this
subsection.

######  Theorem 4.7.

Assume that the objective functions $f_{i}$ are convex with gradients $\nabla
f_{i}$ that are Lipschitz continuous on bounded sets. Then, solution
$x:[0,+\infty)\to\mathcal{H}$ of (13) with arbitrary initial conditions
$x^{0},v^{0}\in\mathcal{H}$ that remains bounded converges weakly to a Pareto
critical point of (8).

###### Proof.

For $z\in S$ define

$\displaystyle h_{z}(t)\coloneqq\frac{1}{2}\lVert x(t)-z\rVert^{2}.$

Using the chain rule, we compute the first and the second derivative of
$h_{z}(t)$ as

$\displaystyle\dot{h}_{z}(t)=\langle
x(t)-z,\dot{x}(t)\rangle,\quad\ddot{h}_{z}(t)=\langle
x(t)-z,\ddot{x}(t)\rangle+\lVert\dot{x}(t)\rVert^{2}.$

For a fixed $t\in(0,+\infty)$, write

$\displaystyle\alpha\dot{h}_{z}(t)+\ddot{h}_{z}(t)=\langle\ddot{x}(t)+\alpha\dot{x}(t),x(t)-z\rangle+\lVert\dot{x}(t)\rVert^{2}.$

Using the definition of (IMOG’), we can write
$\ddot{x}(t)+\alpha\dot{x}(t)=-\sum_{i=1}^{m}\theta_{i}\nabla f_{i}(x(t))$ for
some weights $\theta\in\Delta_{m}$. Then, we write

$\displaystyle\alpha\dot{h}_{z}(t)+\ddot{h}_{z}(t)=\sum_{i=1}^{m}\theta_{i}\langle\nabla
f_{i}(x(t)),z-x(t)\rangle+\lVert\dot{x}(t)\rVert^{2}.$ (24)

Proposition 4.1 gives for all $i=1,\dots,m$

$\displaystyle\mathcal{E}_{i}(t)=f_{i}(x(t))+\frac{1}{2}\lVert\dot{x}(t)\rVert^{2}\geq\mathcal{E}_{i}^{\infty}\geq
f_{i}(z)\geq f_{i}(x(t))+\langle\nabla f_{i}(x(t)),z-x(t)\rangle,$

and therefore

$\displaystyle\sum_{i=1}^{m}\theta_{i}\langle\nabla
f_{i}(x(t)),z-x(t)\rangle\leq\frac{1}{2}\lVert\dot{x}(t)\rVert^{2}.$ (25)

Combining inequalities (24) and (25) we get

$\displaystyle\alpha\dot{h}_{z}(t)+\ddot{h}_{z}(t)\leq\frac{3}{2}\lVert\dot{x}(t)\rVert^{2}.$

By Proposition 4.2, we know $\lVert\dot{x}(\cdot)\rVert^{2}\in
L^{1}([0,+\infty))$. Then, Lemma 4.6 states that $\lim_{t\to+\infty}h_{z}(t)$
exists. In addition, we know that every weak sequential cluster point of
$x(t)$ belongs to $S$ by the weak lower semicontinuity of the objective
functions $f_{i}$. Then, we can use Opial’s Lemma 5.5 to prove that $x(t)$
converges weakly to an element in $S$. Let $x^{\infty}$ be the weak limit of
$x(t)$. Then, by Proposition 4.2, there exists a monotonically increasing
unbounded sequence $(t_{k})_{k\geq 0}$ with
$\operatorname*{proj}_{C(x(t_{k}))}(0)\to 0$ for $k\to+\infty$. Since
$x(t_{k})$ converges weakly to $x^{\infty}$, Lemma 4.3 states that
$x^{\infty}$ is Pareto critical. ∎

##  5 An Inertial Multiobjective Gradient Algorithm

In this section, we derive an inertial first order method for multiobjective
optimization problems from an explicit discretization of the differential
equation (IMOG’). We write the system (IMOG’) in the equivalent form

$\displaystyle\alpha\dot{x}(t)+\operatorname*{proj}_{C(x(t))+\ddot{x}(t)}(0),$

and use the following discretization of the differential equation

$\displaystyle\begin{split}\alpha\frac{x^{k+1}-x^{k}}{h}+\operatorname*{proj}_{C(x^{k})+\frac{x^{k+1}-2x^{k}+x^{k-1}}{h^{2}}}(0)=0,\\\
\alpha
h(x^{k+1}-x^{k})+\operatorname*{proj}_{h^{2}C(x^{k})+x^{k+1}-2x^{k}+x^{k-1}}(0)=0.\end{split}$

Lemma A.2 states that $x^{k+1}$ is uniquely defined as

$\displaystyle\begin{split}x^{k+1}&=-\left(\frac{1}{1+\alpha
h}\operatorname*{proj}_{h^{2}C(x^{k})-2x^{k}+x^{k-1}}(-x^{k})-\frac{\alpha
h}{1+\alpha h}x^{k}\right)\\\ &=-\left(\frac{1}{1+\alpha
h}\left[-x^{k}+\operatorname*{proj}_{h^{2}C(x^{k})-x^{k}+x^{k-1}}(0)\right]-\frac{\alpha
h}{1+\alpha h}x^{k}\right)\\\ &=x^{k}-\frac{1}{1+\alpha
h}\operatorname*{proj}_{h^{2}C(x^{k})-x^{k}+x^{k-1}}(0).\end{split}$ (26)

Therefore, $x^{k+1}$ can be written as

$\displaystyle x^{k+1}=x^{k}+\frac{1}{1+\alpha
h}(x^{k}-x^{k-1})-\frac{h^{2}}{1+\alpha h}\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k}),$ (27)

where $\theta^{k}\in\Delta_{m}$ is the solution to the quadratic optimization
problem

$\displaystyle\begin{split}\min_{\theta\in\mathbb{R}^{m}}\,\left\lVert
h^{2}\left(\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(x^{k})\right)-(x^{k}-x^{k-1})\right\rVert^{2}\,\,\operatorname{s.t.}\,\theta\geq
0\,\,\text{and}\,\,\sum_{i=1}^{m}\theta_{i}=1.\end{split}$ (28)

In the following subsection, we analyze the asymptotic behavior of the
sequence $(x^{k})_{k\geq 0}$ that is defined by equations (27) and (28).

###  5.1 Asymptotic Analysis

The asymptotic analysis of the sequence $(x^{k})_{k\geq 0}$ defined by (27)
and (28) works surprisingly similar to the asymptotic analysis of the
trajectories $x(t)$ of the differential equation (IMOG’). We start by proving
that the sequence $(x^{k})_{k\geq 0}$ satisfies a dissipative property. To
this end, we introduce the following preparatory lemma.

######  Lemma 5.1.

Let $(x^{k})_{k\geq 0}$ be defined by (27) and (28) with
$x^{0}=x^{1}\in\mathcal{H}$ and $\alpha,h>0$. Then, for all $i=1,\dots,m$ it
holds that

$\displaystyle\langle\nabla
f_{i}(x^{k}),x^{k+1}-x^{k}\rangle\leq-\frac{\alpha}{h}\lVert
x^{k+1}-x^{k}\rVert^{2}+\frac{1}{2h^{2}}\left[\lVert
x^{k}-x^{k-1}\rVert^{2}-\lVert x^{k+1}-x^{k}\rVert^{2}\right].$

###### Proof.

Using the variational characterization of the convex projection in the
identity (5), we get for all $i=1,\dots,m$,

$\displaystyle\langle\alpha h(x^{k+1}-x^{k})+h^{2}\nabla
f_{i}(x^{k})+(x^{k+1}-x^{k})-(x^{k}-x^{k-1}),\alpha
h(x^{k+1}-x^{k})\rangle\leq 0,$

which can be rearranged into

$\displaystyle\begin{split}&\langle\nabla f_{i}(x^{k}),x^{k+1}-x^{k}\rangle\\\
\leq&-\left(\frac{1}{h^{2}}+\frac{\alpha}{h}\right)\lVert
x^{k+1}-x^{k}\rVert+\frac{1}{h^{2}}\langle
x^{k+1}-x^{k},x^{k}-x^{k-1}\rangle\\\ \leq&-\frac{\alpha}{h}\lVert
x^{k+1}-x^{k}\rVert^{2}+\frac{1}{2h^{2}}\left[\lVert
x^{k}-x^{k-1}\rVert^{2}-\lVert x^{k+1}-x^{k}\rVert^{2}\right].\end{split}$
(29)

∎

Using Lemma 5.1, we show that there exists an energy sequence which can be
seen as a discretization of the energy function defined in Proposition 4.2.

######  Proposition 5.2.

Assume that the gradients $\nabla f_{i}$ of the objective functions are
globally $L$-Lipschitz continuous for all $i=1,\dots,m$ and further assume
$Lh<2\alpha$. Then

$\displaystyle\mathcal{E}_{i,k}\coloneqq f_{i}(x^{k})+\frac{1}{2h^{2}}\lVert
x^{k}-x^{k-1}\rVert^{2},$ (30)

is monotonically decreasing.

###### Proof.

We start with investigating the difference

$\displaystyle\mathcal{E}_{i,k+1}-\mathcal{E}_{i,k}=f_{i}(x^{k+1})-f_{i}(x^{k})+\frac{1}{2h^{2}}\left[\lVert
x^{k+1}-x^{k}\rVert-\lVert x^{k}-x^{k-1}\rVert\right].$

Using that $f_{i}$ is convex with $L$-Lipschitz continuous gradient, we
estimate the expression above by

$\displaystyle\leq\langle\nabla
f_{i}(x^{k}),x^{k+1}-x^{k}\rangle+\frac{L}{2}\lVert
x^{k+1}-x^{k}\rVert^{2}+\frac{1}{2h^{2}}\left[\lVert
x^{k+1}-x^{k}\rVert-\lVert x^{k}-x^{k-1}\rVert\right].$

Using Lemma 5.1, we estimate this term by

$\displaystyle\leq$
$\displaystyle\left(\frac{L}{2}-\frac{\alpha}{h}\right)\lVert
x^{k+1}-x^{k}\rVert^{2}-\frac{1}{2h^{2}}\lVert
x^{k+1}-2x^{k}+x^{k-1}\rVert^{2}.$

For $hL<2\alpha$ it holds that $\left(\frac{L}{2}-\frac{\alpha}{h}\right)<0$
and we get

$\displaystyle\mathcal{E}_{i,k+1}-\mathcal{E}_{i,k}\leq\left(\frac{L}{2}-\frac{\alpha}{h}\right)\lVert
x^{k+1}-x^{k}\rVert^{2}-\frac{1}{2h^{2}}\lVert
x^{k+1}-2x^{k}+x^{k-1}\rVert^{2},$ (31)

which completes the proof. ∎

The following corollary is an immediate consequence of Proposition 5.2.

######  Corollary 5.3.

Assume all conditions of Proposition 5.2 are met. Then, for all $i=1,\dots,m$
and all $k\geq 1$ it holds that $f_{i}(x^{k})\leq f_{i}(x^{0})$.

Corollary 5.3 hints at a condition that guarantees that the sequence
$(x^{k})_{k\geq 0}$ remains bounded. If the level set
$\mathcal{L}_{i}(f_{i}(x^{0}))\coloneqq\\{x\in\mathcal{H}\,:\,f_{i}(x)\leq
f_{i}(x^{0})\\}$ of one objective function $f_{i}$ is bounded, the sequence
$(x^{k})_{k\geq 0}$ remains bounded. In the following proposition we collect
some immediate consequences of Proposition 5.2.

######  Proposition 5.4.

Assume the gradients $\nabla f_{i}$ of the objective functions are
$L$-Lipschitz continuous on a bounded set containing the sequence
$(x^{k})_{k\geq 0}$, that is defined by equations (27) and (28). Assume
$Lh<2\alpha$, then for all $i=1,\dots,m$ the following statements hold.

  1. i)

$\mathcal{E}_{i,k}\to\mathcal{E}_{i}^{\infty}>-\infty$ as $k\to+\infty$

  2. ii)

$\sum_{k=0}^{\infty}\lVert x^{k+1}-x^{k}\rVert^{2}<+\infty$

  3. iii)

$f_{i}(x^{k})\to\mathcal{E}_{i}^{\infty}$ as $k\to+\infty$

###### Proof.

i) Proposition 5.2 states that $\mathcal{E}_{i,k}$ is monotonically
decreasing. Therefore, $\mathcal{E}_{i,k}\to\mathcal{E}_{i}^{\infty}$ holds.
We have to show that $\mathcal{E}_{i}^{\infty}>-\infty$. Since the objective
functions $f_{i}$ have Lipschitz continuous gradients on a bounded set
containing $(x^{k})_{k\geq 0}$, it follows by the mean value theorem that
$f_{i}$ is bounded on this sets and in particular on $(x^{k})_{k\geq 0}$.
Therefore, we conclude

$\displaystyle\mathcal{E}_{i}^{\infty}=\lim_{k\to+\infty}f_{i}(x^{k})+\frac{1}{2h^{2}}\lVert
x^{k}-x^{k-1}\rVert^{2}\geq\liminf_{k\to+\infty}f_{i}(x^{k})>-\infty.$

ii) From the inequality (31) we immediately follow

$\displaystyle\begin{split}&\mathcal{E}_{i,K+1}-\mathcal{E}_{i,1}=\sum_{k=1}^{K}\left(\mathcal{E}_{i,k+1}-\mathcal{E}_{i,k}\right)\\\
\leq&\sum_{k=1}^{K}\left(\frac{L}{2}-\frac{\alpha}{h}\right)\lVert
x^{k+1}-x^{k}\rVert^{2}-\frac{1}{h^{2}}\sum_{k=1}^{K}\lVert
x^{k+1}-2x^{k}+x^{k-1}\rVert^{2}.\end{split}$ (32)

Since $Lh<2\alpha$, it holds that
$\left(\frac{\alpha}{h}-\frac{L}{2}\right)>0$ and therefore we get for all
$K\geq 1$

$\displaystyle\left(\frac{\alpha}{h}-\frac{L}{2}\right)\sum_{k=1}^{K}\lVert
x^{k+1}-x^{k}\rVert^{2}\leq\mathcal{E}_{i,1}-\mathcal{E}_{i,K+1}.$

From part i), we know that the right hand side converges which completes the
proof of ii).

iii) Since $\mathcal{E}_{i,k}\to\mathcal{E}_{i}^{\infty}$ and $\lVert
x^{k+1}-x^{k}\rVert^{2}\to 0$, it follows that
$f_{i}(x^{k})\to\mathcal{E}_{i}^{\infty}$. ∎

We use the following discrete version of Opial’s Lemma to prove that
$(x^{k})_{k\geq 0}$ converges weakly to a Pareto critical point of (8).

######  Lemma 5.5 (Opial’s Lemma [28]).

Let $S\subset\mathcal{H}$ be nonempty and let $(x_{k})_{k\geq 0}$ be a
sequence in $\mathcal{H}$ that satisfies the following conditions.

  1. i)

For all $z\in S$ $\lim_{k\to+\infty}\lVert x^{k}-z\rVert$ exists.

  2. ii)

Every weak sequential cluster point of $(x_{k})_{k\geq 0}$ belongs to $S$.

Then, it follows that $(x_{k})_{k\geq 0}$ converges weakly to an element in
$S$.

We will use Opial’s Lemma on the set

$\displaystyle
S\coloneqq\left\\{z\in\mathcal{H}\,:\,f_{i}(z)\leq\lim_{k\to+\infty}f_{i}(x_{k})\,\,\text{for
all }i=1,\dots,m,\right\\}.$ (33)

######  Theorem 5.6.

Assume the gradients $\nabla f_{i}$ of the objective functions are
$L$-Lipschitz continuous on a bounded set containing the sequence
$(x^{k})_{k\geq 0}$, defined by (27) and (28) and further assume $Lh<2\alpha$.
Then, $(x^{k})_{k\geq 0}$ converges weakly to a Pareto critical point of (8).

###### Proof.

We show that $(x^{k})_{k\geq 0}$ satisfies Opial’s Lemma for the set $S$
defined by (33). We start by showing a quasi Fejér property of the sequence
$(x_{k})_{k\geq 0}$. For a fixed $z\in S$, define the sequence

$\displaystyle h_{k}\coloneqq\frac{1}{2}\lVert x^{k}-z\rVert^{2}.$

It is easy to check that

$\displaystyle h_{k+1}$ $\displaystyle=h_{k}+\langle
x^{k+1}-x^{k},x^{k}-z\rangle+\frac{1}{2}\lVert x^{k+1}-x^{k}\rVert^{2}.$

Proposition 5.4 guarantees the monotonicity of $\mathcal{E}_{i,k}$. Since
$z\in S$, from the convexity of $f_{i}$ we can deduce for all $i=1,\dots,m$
that

$\displaystyle\mathcal{E}_{i,k}=f_{i}(x^{k})+\frac{1}{2h^{2}}\lVert
x^{k}-x^{k-1}\rVert^{2}\geq\mathcal{E}_{i}^{\infty}\geq f_{i}(z)\geq
f_{i}(x^{k})+\langle\nabla f_{i}(x^{k}),z-x^{k}\rangle,$

and therefore

$\displaystyle\langle\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(x^{k}),z-x^{k}\rangle\leq\frac{1}{2h^{2}}\lVert
x^{k}-x^{k-1}\rVert^{2}.$

Using this inequality we can show

$\displaystyle\frac{h^{2}}{1+\alpha
h}\langle\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k}),z-x^{k}\rangle=\langle x^{k}-x^{k+1}-\frac{1}{1+\alpha
h}(x^{k}-x^{k-1}),z-x^{k}\rangle$ $\displaystyle=$ $\displaystyle\langle
x^{k+1}-x^{k},x^{k}-z\rangle-\frac{1}{1+\alpha h}\langle
x^{k-1}-x^{k},x^{k}-z\rangle\leq\frac{1}{2(1+\alpha h)}\lVert
x^{k}-x^{k-1}\rVert^{2},$

which leads to the inequality

$\displaystyle\langle x^{k+1}-x^{k},x^{k}-z\rangle\leq\frac{1}{1+\alpha
h}\langle x^{k-1}-x^{k},x^{k}-z\rangle+\frac{1}{2(1+\alpha h)}\lVert
x^{k}-x^{k-1}\rVert^{2}.$

We use this inequality to show

$\displaystyle h_{k+1}-h_{k}=$ $\displaystyle\langle
x^{k+1}-x^{k},x^{k}-z\rangle+\frac{1}{2}\lVert x^{k+1}-x^{k}\rVert^{2}$
$\displaystyle\leq$ $\displaystyle\frac{1}{1+\alpha h}\langle
x^{k-1}-x^{k},x^{k}-z\rangle+\frac{1}{2}\lVert
x^{k+1}-x^{k}\rVert^{2}+\frac{1}{2(1+\alpha h)}\lVert x^{k}-x^{k-1}\rVert^{2}$
$\displaystyle=$ $\displaystyle\frac{1}{1+\alpha
h}\left[h_{k}-h_{k-1}+\frac{1}{2}\lVert
x^{k}-x^{k-1}\rVert^{2}\right]+\frac{1}{2}\lVert x^{k+1}-x^{k}\rVert^{2}$
$\displaystyle+\frac{1}{2(1+\alpha h)}\lVert x^{k}-x^{k-1}\rVert^{2}$
$\displaystyle\leq$ $\displaystyle\frac{1}{1+\alpha
h}(h_{k}-h_{k-1})+\frac{1}{2}\lVert x^{k+1}-x^{k}\rVert^{2}+\frac{1}{1+\alpha
h}\lVert x^{k}-x^{k-1}\rVert^{2}.$

Defining $\theta_{k}\coloneqq h_{k+1}-h_{k}$,
$\delta_{k}\coloneqq\frac{1}{1+\alpha h}\lVert
x^{k}-x^{k-1}\rVert^{2}+\frac{1}{2}\lVert x^{k+1}-x^{k}\rVert^{2}$ and
$a\coloneqq\frac{1}{1+\alpha h}$, we can therefore conclude

$\displaystyle\theta_{k+1}\leq a\theta_{k}+\delta_{k}.$

Proposition 5.4 states that $\sum_{k=1}^{\infty}\delta_{k}<+\infty$.
Therefore, we can use Theorem 2.1 in [29] or Theorem 3.1 in [4] to show that
$h_{k}$ converges. To use Opial’s Lemma, we also have to show that all weak
sequential cluster points of $(x_{k})_{k\geq 0}$ belong to $S$. Since the
sequence $(x_{k})_{k\geq 0}$ is bounded, it possesses at least one sequential
cluster point that we denote by $x^{\infty}$ and a subsequence
$(x_{k_{l}})_{l\geq 0}$ that converges weakly to $x^{\infty}$. Since $f_{i}$
is convex and continuous, it is also weakly lower semicontinuous and it
follows that for all $i=1,\dots,m$

$\displaystyle
f_{i}(x^{\infty})\leq\liminf_{l\to+\infty}f_{i}(x_{k_{l}})=\lim_{k\to+\infty}f_{i}(x_{k}),$

where the equality follows from the fact that the limit exists. Therefore,
$x^{\infty}\in S$ and hence $S$ is nonempty. Then, Opial’s Lemma 5.5 states
that $(x^{k})_{k\geq 0}$ converges weakly to an element in $S$ that we denote
by $x^{\infty}$. We will show that each weak sequential cluster point of
$(x^{k})_{k\geq 0}$ is Pareto critical. By the definition of the sequence
$(x^{k})_{k\geq 0}$ in (27), it holds that

$\displaystyle\sum_{k=1}^{\infty}\left\lVert\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k})\right\rVert^{2}$
$\displaystyle=\sum_{k=1}^{\infty}\left\lVert\frac{1+\alpha
h}{h^{2}}(x^{k+1}-x^{k})+\frac{1}{h^{2}}(x^{k}-x^{k-1})\right\rVert^{2}.$

This sum is finite by part ii) of Proposition 5.4. Thus, we know that the
sequence $g^{k}\coloneqq\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k})\in\operatorname{co}(\nabla f_{i}(x^{k}))$ converges strongly to
zero. Since $x^{k}$ converges weakly to $x^{\infty}$, Lemma 4.3 states that
$0\in\operatorname{co}(\nabla f_{i}(x^{\infty}))$ and hence $x^{\infty}$ is
Pareto critical. ∎

##  6 An Accelerated Multiobjective Gradient Method

In this section, we define a multiobjective gradient method with Nesterov
acceleration based on the inertial method we discussed in the previous
subsection.

###  6.1 The single objective Case

In this subsection we present Nesterov’s method in the single objective
setting and point out its relation to an intertial gradient-like dynamical
system with asymptotically vanishing damping. Consider the problem

$\displaystyle\min_{x\in\mathcal{H}}f(x),$

where $f:\mathcal{H}\to\mathbb{R}$ is convex and differentiable with
$L$-Lipschitz continuous gradient $\nabla f(x)$. For $\alpha\geq 3$,
$0<s\leq\frac{1}{L}$ and $x^{0},x^{1}\in\mathcal{H}$, define the sequence
$(x^{k})_{k\geq 0}$ by

$\displaystyle\left.\begin{array}[]{rl}y^{k}&=x^{k}+\frac{k-1}{k+\alpha-1}(x^{k}-x^{k-1}),\\\
x^{k+1}&=y^{k}-s\nabla f(y^{k})\end{array}\right\\}\,\,\text{for}\,\,k\geq 1.$
(36)

If $\operatorname*{arg\,min}f\neq\emptyset$ it can be shown that
$f(x^{k})-\min_{x\in\mathcal{H}}f(x)=\mathcal{O}(k^{-2})$ and that $\lVert
x^{k+1}-x^{k}\rVert=\mathcal{O}(k^{-1})$. For $\alpha>3$, it holds that
$f(x^{k})-\min_{x\in\mathcal{H}}f(x)=o(k^{-2})$, $\lVert
x^{k+1}-x^{k}\rVert=o(k^{-1})$ and that $(x^{k})_{k\geq 0}$ converges weakly
to an element in $\operatorname*{arg\,min}f$ [6]. Nesterov’s method is related
to the following gradient-like dynamical system with asymptotically vanishing
damping

$\displaystyle\ddot{x}(t)+\frac{\alpha}{t}\dot{x}(t)+\nabla f(x(t))=0.$ (37)

The algorithm (36) can be derived as a discretization of (37). This relation
is further investigated in [1, 30].

###  6.2 Introducing Nesterov Acceleration in (IMOG’)

We formally define the following gradient-like system with asymptotically
vanishing damping for multiobjective optimization.

$\displaystyle\ddot{x}(t)+\frac{3}{t}\dot{x}(t)+\operatorname*{proj}_{C(x(t))}{(-\ddot{x}(t))}=0.$
(38)

We will not analyze this system in detail but leave this for future work. We
show that an implicit discretization of this system leads to an accelerated
multiobjective gradient method with an improvement convergence rate of the
function values. We equivalently write (38) as

$\displaystyle\frac{3}{t}\dot{x}(t)+\operatorname*{proj}_{C(x(t))+\ddot{x}(t)}(0)=0.$
(39)

Using the same Ansatz as in Section 2 of [1], we show that we can derive the
differential equation (38) from the scheme

$\displaystyle\frac{3}{k}(x^{k+1}-x^{k})+\operatorname*{proj}_{sC(y^{k})+(x^{k+1}-2x^{k}+x^{k-1})}(0)=0,$
(40)

with $y^{k}=x^{k}+\frac{k-1}{k+2}(x^{k}-x^{k-1})$. We divide (40) by
$\sqrt{s}$ to get

$\displaystyle\frac{3}{k}\frac{x^{k+1}-x^{k}}{\sqrt{s}}+\operatorname*{proj}_{\sqrt{s}C(y^{k})+\frac{x^{k+1}-2x^{k}+x^{k-1}}{\sqrt{s}}}(0)=0.$
(41)

We use the Ansatz $x_{k}\approx x(k\sqrt{s})$ for some smooth curve $x(t)$
defined for all $t\geq 0$. Write $k=\frac{t}{\sqrt{s}}$. When the step size
$s$ goes to zero $X(t)\approx x_{\frac{t}{\sqrt{s}}}=x_{k}$ and $X(t)\approx
x_{\frac{t+\sqrt{s}}{\sqrt{s}}}=x_{k+1}$. Then, Taylor expansion gives

$\displaystyle\frac{x^{k+1}-x^{k}}{\sqrt{s}}=\dot{x}(t)+\frac{1}{2}\ddot{x}(t)\sqrt{s}+o(\sqrt{s}),\quad\frac{x^{k}-x^{k-1}}{\sqrt{s}}=\dot{x}(t)-\frac{1}{2}\ddot{x}(t)\sqrt{s}+o(\sqrt{s}),$
(42)

and hence

$\displaystyle\frac{x^{k+1}-2x^{k}+x^{k-1}}{\sqrt{s}}=\ddot{x}(t)\sqrt{s}+o(\sqrt{s}).$
(43)

For all $i=1,\dots,m$, we have $\sqrt{s}\nabla f_{i}(y^{k})=\sqrt{s}\nabla
f_{i}(x(t))+o(\sqrt{s})$. Since the convex projection depends in a well-
behaved manner on the convex set we project onto, we get

$\displaystyle\operatorname*{proj}_{\sqrt{s}C(y^{k})+\frac{x^{k+1}-2x^{k}+x^{k-1}}{\sqrt{s}}}(0)=\sqrt{s}\operatorname*{proj}_{C(x(t))+\ddot{x}(t)}(0)+o(\sqrt{s}).$
(44)

Combining (42), (43) and (44), we get from (41)

$\displaystyle\frac{3\sqrt{s}}{t}\left(\dot{x}(t)+\frac{1}{2}\ddot{x}(t)\sqrt{s}+o(\sqrt{s})\right)+\sqrt{s}\operatorname*{proj}_{C(x(t))+\ddot{x}(t)}(0)+o(\sqrt{s})=0.$

Comparing the coefficients of $\sqrt{s}$, we obtain

$\displaystyle\frac{3}{t}\dot{x}(t)+\operatorname*{proj}_{C(x(t))+\ddot{x}(t)}(0)=0.$

We have shown that the differential equation (38) can be derived from the
scheme (40). Using Lemma A.1 on (40), we get that $x^{k+1}$ is uniquely
defined as

$\displaystyle x^{k+1}$
$\displaystyle=-\left(\frac{k}{k+3}\operatorname*{proj}_{sC(y^{k})-2x^{k}+x^{k-1}}(-x^{k})-\frac{3}{k+3}x^{k}\right)$
$\displaystyle=x^{k}-\frac{k}{k+3}\operatorname*{proj}_{sC(y^{k})-x^{k}-x^{k-1}}(0).$

The last term can be witten as

$\displaystyle
x^{k}+\frac{k}{k+3}(x^{k}-x^{k-1})-\frac{sk}{k+3}\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(y^{k}),$

where $\theta^{k}\in\mathbb{R}^{m}$ is a solution to the quadratic
optimization problem

$\displaystyle\begin{split}\min_{\theta\in\mathbb{R}^{m}}\,\left\lVert
s\left(\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k})\right)-(x^{k}-x^{k-1})\right\rVert^{2}\,\,\operatorname{s.t.}\,\theta\geq
0\,\,\text{and}\,\,\sum_{i=1}^{m}\theta_{i}=1.\end{split}$ (45)

We want to drop the factor $\frac{k}{k+3}$ in front of the term
$\sum_{i=1}^{m}\theta_{i}^{k}\nabla f_{i}(y^{k})$ to get a method that more
closely resembles (36). In addition, we perform a shift of the index $k$ to
transform $\frac{k}{k+3}$ into $\frac{k-1}{k+2}$. The final method we define
in this subsection can then be defined as follows. Let
$x^{0}=x^{1}\in\mathcal{H}$ and $s>0$. Define the scheme

$\displaystyle\begin{split}\left.\begin{array}[]{c}y^{k}=x^{k}+\frac{k-1}{k+2}(x^{k}-x^{k-1}),\\\
x^{k+1}=y^{k}-s\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(y^{k})\end{array}\right\\}\,\,\text{ for }\,\,k\geq 1,\end{split}$ (46)

where in each step $\theta^{k}\in\mathbb{R}^{m}$ is a solution to the
quadratic optimization problem

$\displaystyle\begin{split}\min_{\theta\in\mathbb{R}^{m}}\,\left\lVert
s\left(\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k})\right)-\frac{k-1}{k+2}(x^{k}-x^{k-1})\right\rVert^{2}\,\,\operatorname{s.t.}\,\theta\geq
0\,\,\text{and}\,\,\sum_{i=1}^{m}\theta_{i}=1.\end{split}$ (47)

The fact that we have to transform the quadratic optimization problem from
(45) into (47) is an observation from the proof of Proposition 6.1. The
presented method is still asymptotically equivalent to the scheme defined by
(40). We summarize the defined method in Algorithm 1 for later references.

Algorithm 1 Accelerated multiobjective gradient method

1:Choose $x^{0}=x^{1}\in\mathcal{H}$, $0<s\leq\frac{1}{L}$ and set $k=1$.

2:Set $y^{k}=x^{k}+\frac{k-1}{k+2}(x^{k}-x^{k-1})$.

3:Compute $\theta^{k}\in\mathbb{R}^{m}$ by solving
$\displaystyle\begin{split}\min_{\theta\in\mathbb{R}^{m}}\,\left\lVert
s\left(\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k})\right)-\frac{k-1}{k+2}(x^{k}-x^{k-1})\right\rVert^{2}\,\,\operatorname{s.t.}\,\theta\geq
0\,\,\text{and}\,\,\sum_{i=1}^{m}\theta_{i}=1.\end{split}$

4:Set $x^{k+1}=y^{k}-s\sum_{i=1}^{m}\theta_{i}^{k}\nabla f_{i}(y^{k})$

5:if stopping condition is true then

6: Stop.

7:else

8: Update $k\leftarrow k+1$ and go to step 1.

9:end if

###  6.3 A Dissipative Property

We start our investigations of Algorithm 1 with an energy estimate analogous
to Proposition 5.2 for the inertial method.

######  Proposition 6.1.

Assume that the gradients $\nabla f_{i}$ of the objective functions are
globally $L$-Lipschitz continuous for all $i=1,\dots,m$ and further assume
$sL\leq 1$. Define for all $k\geq 1$ the energy sequence

$\displaystyle\mathcal{E}_{i,k}\coloneqq f_{i}(x^{k})+\frac{1}{2s}\lVert
x^{k}-x^{k-1}\rVert^{2}.$

For all $k\geq 1$, it holds that

$\displaystyle\mathcal{E}_{i,k+1}-\mathcal{E}_{i,k}\leq-\frac{1}{2s}\frac{3}{k+2}\lVert
x^{k}-x^{k-1}\rVert^{2}.$

###### Proof.

From the definition of $x^{k}$ and $y^{k}$ in (46) we get

$\displaystyle
x^{k+1}-x^{k}+\operatorname*{proj}_{sC(y^{k})-\frac{k-1}{k+2}(x^{k}-x^{k-1})}(0)$
$\displaystyle=0.$

Hence, for all $i=1,\dots,m$ it holds that

$\displaystyle\langle x^{k+1}-x^{k}+s\nabla
f_{i}(y^{k})-\frac{k-1}{k+2}(x^{k}-x^{k-1}),x^{k+1}-x^{k}\rangle\leq 0,$

from which we follow

$\displaystyle s\langle\nabla f_{i}(y^{k}),x^{k+1}-x^{k}\rangle\leq-\lVert
x^{k+1}-x^{k}\rVert^{2}+\frac{k-1}{k+2}\langle
x^{k+1}-x^{k},x^{k}-x^{k-1}\rangle$ $\displaystyle=$
$\displaystyle-\frac{3}{k+2}\lVert
x^{k+1}-x^{k}\rVert-\frac{1}{2}\frac{k-1}{k+2}\lVert
x^{k+1}-2x^{k}+x^{k-1}\rVert^{2}$
$\displaystyle+\frac{1}{2}\frac{k-1}{k+2}\left[\lVert
x^{k}-x^{k-1}\rVert-\lVert x^{k+1}-x^{k}\rVert\right].$

Writing out the definition of $y^{k}$, one can easily verify that

$\displaystyle\lVert x^{k+1}-y^{k}\rVert^{2}\leq\frac{k-1}{k+2}\lVert
x^{k+1}-2x^{k}+x^{k-1}\rVert^{2}+\frac{3}{k+2}\lVert x^{k+1}-x^{k}\rVert^{2}.$

Combining the inequalities above and using $sL\leq 1$ we get

$\displaystyle s(f_{i}(x^{k+1})-f_{i}(x^{k}))\leq s\langle\nabla
f_{i}(y^{k}),x^{k+1}-x^{k}\rangle+\frac{1}{2}\lVert x^{k+1}-y^{k}\rVert^{2}$
$\displaystyle\leq$ $\displaystyle-\frac{1}{2}\frac{3}{k+2}\lVert
x^{k+1}-x^{k}\rVert^{2}+\frac{1}{2}\frac{k-1}{k+2}\left[\lVert
x^{k}-x^{k-1}\rVert-\lVert x^{k+1}-x^{k}\rVert\right]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\lVert x^{k}-x^{k-1}\rVert^{2}-\lVert
x^{k+1}-x^{k}\rVert^{2}\right]-\frac{1}{2}\frac{3}{k+2}\lVert
x^{k}-x^{k-1}\rVert^{2},$

which completes the proof. ∎

######  Corollary 6.2.

Let $(x^{k})_{k\geq 0}$ be a sequence defined by (46). Then, it holds that for
all $k\geq 0$ and all $i=1,\dots,m$

$\displaystyle f_{i}(x^{k})\leq f_{i}(x^{0}).$

###  6.4 Convergence of Function Values with Rate $\mathcal{O}(k^{-2})$

The proof in this section relies on the proof by Fukuda, Tanabe and Yamashita
[9] for their accelerated gradient method and the proof of Attouch and
Peypouquet [6] for the single objective case. The following definition is
aligned with [22] and the concept of merit functions that was introduced in
[22] and further utilized in [31, 9]. For $z\in\mathcal{H}$ define

$\displaystyle\sigma_{k}(z)\coloneqq\min_{i=1,\dots,m}f_{i}(x^{k})-f_{i}(z).$

######  Lemma 6.3.

It holds that

$\displaystyle\sigma_{k+1}(z)$ $\displaystyle\leq-\frac{1}{s}\langle
x^{k+1}-y^{k},y^{k}-z\rangle-\frac{1}{2s}\lVert x^{k+1}-y^{k}\rVert^{2}.$

###### Proof.

The objective functions $f_{i}$ are convex with $L$-Lipschitz continuous
gradients. Therefore, for all $i=1,\dots,m$ it holds that

$\displaystyle\begin{split}&f_{i}(x^{k+1})-f_{i}(z)\leq
f_{i}(x^{k+1})-f_{i}(y^{k})+f_{i}(y^{k})-f_{i}(z)\\\ \leq&\langle\nabla
f_{i}(y^{k}),x^{k+1}-y^{k}\rangle+\frac{L}{2}\lVert
x^{k+1}-y^{k}\rVert^{2}+\langle\nabla f_{i}(y^{k}),y^{k}-z\rangle.\end{split}$
(48)

The definition of $\sigma_{k}(z)$ gives

$\displaystyle\sigma_{k+1}(z)=\min_{i=1,\dots,m}f_{i}(x^{k+1})-f_{i}(z)\leq\sum_{i=1}^{m}\theta_{i}^{k}\left(f_{i}(x^{k+1})-f_{i}(z)\right).$
(49)

Combining (48) and (49) and using $\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(y^{k})=\frac{1}{s}(y^{k}-x^{k+1})$ we get the desired inequality. ∎

We want to find a similar inequality for the expression
$f_{i}(x^{k+1})-f_{i}(x^{k})$. To this end, we introduce the following lemma.

######  Lemma 6.4.

Define the optimization problem

$\displaystyle\begin{split}\min_{(v,\alpha)\in\mathcal{H}\times\mathbb{R}}&\Phi(v,\alpha)\coloneqq\frac{1}{2}\lVert
sv+(y^{k}-x^{k})\rVert^{2}+\alpha\\\
\operatorname{s.t.}\,\,&g_{i}(v,\alpha)\coloneqq\langle s\nabla
f_{i}(y^{k})-(y^{k}-x^{k}),sv+(y^{k}-x^{k})\rangle-\alpha\leq 0.\end{split}$
(50)

Then, it holds that the dual problem to this problem is the quadratic problem
(47). An optimal solution $\theta^{*}$ to (47) satisfies

$\displaystyle\langle s\sum_{i=1}^{m}\theta_{i}^{*}\nabla
f_{i}(y^{k}),x^{k+1}-x^{k}\rangle=\max_{i=1,\dots,m}\langle s\nabla
f_{i}(y^{k}),x^{k+1}-x^{k}\rangle.$

###### Proof.

Since $\mathcal{H}$ is potentially infinite-dimensional, we need duality
statements for infinite-dimensional constrained optimization problems. The
statements we use in this proof can be found in Sections 8.3 to 8.6 of [32].
Since the optimization problem (50) has a fairly simple structure, we will not
recite every result we use. The duality between (50) and (47) follows from a
straightforward computation. Since the objective function $\Phi(v,\alpha)$ of
(50) is convex and all constraints $g_{i}(v,\alpha)$ are linear, strong
duality holds. Hence a KKT point
$((v^{*},\alpha^{*}),\theta^{*})\in(\mathcal{H}\times\mathbb{R})\times\mathbb{R}^{m}$
of problem (50) yields a solution to (47). From the KKT conditions for (50) we
get that

$\displaystyle v^{*}=-s\sum_{i=1}^{m}\theta_{i}^{*}\nabla f_{i}(y^{k}).$

For all $i=1,\dots,m$ it holds that $g_{i}(v,\alpha)\leq 0$ and hence

$\displaystyle\langle s\nabla
f_{i}(y^{k})-(y^{k}-x^{k}),sv+(y^{k}-x^{k})\rangle\leq\alpha.$

By the complementarity of $\theta_{i}^{*}$ and $g_{i}(v^{*},\alpha^{*})$ we
get

$\displaystyle\langle s\sum_{i=1}^{m}\theta_{i}^{*}\nabla
f_{i}(y^{k})-(y^{k}-x^{k}),sv^{*}+(y^{k}-x^{k})\rangle=\alpha^{*}$
$\displaystyle=$ $\displaystyle\max_{i=1,\dots,m}\langle s\nabla
f_{i}(y^{k})-(y^{k}-x^{k}),sv^{*}+(y^{k}-x^{k})\rangle.$

The second equality above follows from the fact that $\theta_{j}^{*}>0$ holds
for at least one $j\in\\{1,\dots,m\\}$ as a consequence of the dual
feasibility.
Using $v^{*}=-\sum_{i=1}^{m}\theta_{i}^{*}\nabla f_{i}(y^{k})$, we get
$sv^{*}=x^{k+1}-y^{k}$ and therefore

$\displaystyle\langle s\sum_{i=1}^{m}\theta_{i}^{*}\nabla
f_{i}(y^{k})-(y^{k}-x^{k}),x^{k+1}-x^{k}\rangle$ $\displaystyle=$
$\displaystyle\max_{i=1,\dots,m}\langle s\nabla
f_{i}(y^{k})-(y^{k}-x^{k}),x^{k+1}-x^{k})\rangle.$

∎

######  Lemma 6.5.

$\displaystyle\sigma_{k+1}(z)-\sigma_{k}(z)$
$\displaystyle\leq-\frac{1}{s}\langle
x^{k+1}-y^{k},y^{k}-x^{k}\rangle-\frac{1}{2s}\lVert x^{k+1}-y^{k}\rVert^{2}.$

###### Proof.

For all $a,b\in\mathbb{R}^{m}$ it holds that

$\displaystyle\left(\min_{i=1,\dots,m}a_{i}\right)-\left(\min_{i=1,\dots,m}b_{i}\right)\leq\max_{i=1,\dots,m}\left(a_{i}-b_{i}\right)$

and therefore for all $z\in\mathcal{H}$

$\displaystyle\sigma_{k+1}(z)-\sigma_{k}(z)\leq\max_{i=1,\dots,m}\left(f_{i}(x^{k+1})-f_{i}(x^{k})\right).$

Using that the objective functions $f_{i}$ are convex with $L$-Lipschitz
continuous gradients and the fact that $sL\leq 1$, we can bound this
expression by

$\displaystyle\leq\max_{i=1,\dots,m}\left(\langle\nabla
f_{i}(y^{k}),x^{k+1}-x^{k}\rangle+\frac{1}{2s}\lVert
x^{k+1}-y^{k}\rVert^{2}\right).$

Now we use Lemma 6.4 and get the equality

$\displaystyle=\sum_{i=1}^{m}\theta_{i}^{k}\langle\nabla
f_{i}(y^{k}),x^{k+1}-x^{k}\rangle+\frac{1}{2s}\lVert x^{k+1}-y^{k}\rVert^{2}.$

From here, we continue by using the definitions of $x^{k}$ and $y^{k}$ from
(46) to get

$\displaystyle=$ $\displaystyle\frac{1}{s}\langle
y^{k}-x^{k+1},x^{k+1}-x^{k}\rangle+\frac{1}{2s}\lVert x^{k+1}-y^{k}\rVert^{2}$
$\displaystyle=$ $\displaystyle-\frac{1}{s}\langle
x^{k+1}-y^{k},y^{k}-x^{k}\rangle-\frac{1}{2s}\lVert x^{k+1}-y^{k}\rVert^{2}.$

∎

######  Theorem 6.6.

The sequence $(x^{k})_{k\geq 0}$ defined by (46) satisfies

$\displaystyle\sigma_{k}(z)\leq\frac{2\left(\lVert x^{1}-z\rVert^{2}+\lVert
x^{2}-z\rVert^{2}\right)}{s(k+1)^{2}}.$

###### Proof.

Lemma 6.3 and Lemma 6.5 state

$\displaystyle\sigma_{k+1}(z)$ $\displaystyle\leq-\frac{1}{s}\langle
x^{k+1}-y^{k},y^{k}-z\rangle-\frac{1}{2s}\lVert
x^{k+1}-y^{k}\rVert^{2}\,\,\text{ and }$
$\displaystyle\sigma_{k+1}(z)-\sigma_{k}(z)$
$\displaystyle\leq-\frac{1}{s}\langle
x^{k+1}-y^{k},y^{k}-x^{k}\rangle-\frac{1}{2s}\lVert x^{k+1}-y^{k}\rVert^{2}.$

Taking a convex combination of the last inequalities with weights
$\frac{2}{k+2}$ and $\frac{k}{k+2}$ yields

$\displaystyle\begin{split}&\sigma_{k+1}(z)-\frac{k}{k+2}\sigma_{k}(z)\\\
\leq&-\frac{1}{s}\langle
x^{k+1}-y^{k},y^{k}-\frac{k}{k+2}x^{k}-\frac{2}{k+2}z\rangle-\frac{1}{2s}\lVert
x^{k+1}-y^{k}\rVert^{2}\\\ =&\frac{1}{s}\langle
x^{k+1}-y^{k},\frac{k}{k+2}(x^{k}-y^{k})+\frac{2}{k+2}(z-y^{k})\rangle-\frac{1}{2s}\lVert
x^{k+1}-y^{k}\rVert^{2}.\end{split}$ (51)

Define

$\displaystyle
z^{k}\coloneqq\frac{k+2}{2}y^{k}-\frac{k}{2}x^{k}=x^{k}+\frac{k-1}{2}(x^{k}-x^{k-1}),$
(52)

and notice that

$\displaystyle\frac{k}{k+2}(y^{k}-x^{k})+\frac{2}{k+2}(y^{k}-z)=\frac{2}{k+2}(z^{k}-z).$
(53)

Using the identity (53) in (51) we get

$\displaystyle\sigma_{k+1}(z)\leq\frac{k}{k+2}\sigma_{k}(z)-\frac{2}{s(k+2)}\langle
x^{k+1}-y^{k},z^{k}-z\rangle-\frac{1}{2s}\lVert x^{k+1}-y^{k}\rVert^{2}.$ (54)

From the definition of $z^{k}$ in (52) one can see that

$\displaystyle z^{k+1}=z^{k}+\frac{k+2}{2}(x^{k+1}-y^{k}).$

Using this identity we can simply compute the squared norm of $\lVert
z^{k+1}-z\rVert^{2}$ as

$\displaystyle\lVert z^{k+1}-z\rVert^{2}=\lVert z^{k}-z\rVert^{2}+(k+2)\langle
z^{k}-z,x^{k+1}-y^{k}\rangle+\left(\frac{k+2}{2}\right)^{2}\lVert
x^{k+1}-y^{k}\rVert^{2}.$

Rearranging this identity and multiplying with $\frac{2}{s(k+2)^{2}}$ yields

$\displaystyle\begin{split}&\frac{2}{s(k+2)^{2}}\left(\lVert
z^{k}-z\rVert^{2}-\lVert z^{k+1}-z\rVert^{2}\right)\\\
=&-\frac{2}{s(k+2)}\langle z^{k}-z,x^{k+1}-y^{k}\rangle-\frac{1}{2s}\lVert
x^{k+1}-y^{k}\rVert^{2}.\end{split}$ (55)

Combining (54) and (55), in total we get

$\displaystyle\sigma_{k+1}(z)\leq\frac{k}{k+2}\sigma_{k}(z)+\frac{4}{2s(k+2)^{2}}\left(\lVert
z^{k}-z\rVert^{2}-\lVert z^{k+1}-z\rVert^{2}\right).$

Multiplying both sides with $(k+2)^{2}$ then yields

$\displaystyle(k+2)^{2}\sigma_{k+1}(z)\leq
k(k+2)\sigma_{k}(z)+\frac{2}{s}\left(\lVert z^{k}-z\rVert^{2}-\lVert
z^{k+1}-z\rVert^{2}\right).$

Using $k(k+2)\leq(k+1)^{2}$ we get

$\displaystyle(k+2)^{2}\sigma_{k+1}(z)-(k+1)^{2}\sigma_{k}(z)\leq\frac{2}{s}\left(\lVert
z^{k}-z\rVert^{2}-\lVert z^{k+1}-z\rVert^{2}\right).$

Summing this inequality from $k=1,\dots,K$, we get for all $z\in\mathcal{H}$

$\displaystyle(K+2)^{2}\sigma_{K+1}(z)\leq\frac{2}{s}\lVert
x^{1}-z\rVert^{2}+4\sigma_{1}(z).$

Similar computations to Lemma 6.3 yield

$\displaystyle\sigma_{1}(z)\leq\frac{1}{2s}\lVert
x^{2}-z\rVert^{2}-\frac{1}{2s}\lVert x^{2}-x^{1}\rVert^{2}.$

Then, for all $k\geq 1$, we obtain

$\displaystyle\sigma_{k}(z)\leq\frac{2\left(\lVert x^{1}-z\rVert^{2}+\lVert
x^{2}-z\rVert^{2}\right)}{s(k+1)^{2}}.$

∎

The theorem above is not straight forward to interpret since we only get
convergence of order $\mathcal{O}(k^{-2})$ for
$\min_{i=1,\dots,m}f_{i}(x^{k})-f_{i}(z)$. This on it’s own does not state
that the vector $f(x^{k})=\left(f_{1}(x^{k}),\dots,f_{m}(x^{k})\right)$
converges to an element of the Pareto front. However we can refine the
statement of Theorem 6.6 in the following way to get a stronger convergence
statement under a weak additional assumption.

######  Theorem 6.7.

Assume in addition to the assumption in Theorem 6.6 that for all
$x\in\mathcal{L}(f(x_{0}))$ there exists an $x^{*}\in\mathcal{L}^{*}\coloneqq
P_{w}\cap\mathcal{L}(f(x_{0}))$ with $f(x^{*})\leq f(x)$ and

$\displaystyle\sup_{f^{*}\in f(\mathcal{L}^{*})}\inf_{x\in
f^{-1}(\\{f^{*}\\})}\lVert x-x^{0}\rVert<+\infty.$ (56)

Then, there exists $R\geq 0$ with

$\displaystyle\sup_{z\in\mathcal{H}}\sigma_{k}(z)\leq\frac{R}{(k+1)^{2}}.$

###### Proof.

Theorem 6.6 gives for all $z\in\mathcal{H}$

$\displaystyle\sigma_{k}(z)\leq\frac{2\left(\lVert x^{1}-z\rVert^{2}+\lVert
x^{2}-z\rVert^{2}\right)}{s(k+1)^{2}}.$

Taking a supremum over this inequality we get

$\displaystyle\sup_{f^{*}\in f(\mathcal{L^{*}})}\inf_{z\in
f^{-1}(f^{*})}\sigma_{k}(z)\leq\sup_{f^{*}\in f(\mathcal{L^{*}})}\inf_{z\in
f^{-1}(f^{*})}\frac{2\left(\lVert x^{1}-z\rVert^{2}+\lVert
x^{2}-z\rVert^{2}\right)}{s(k+1)^{2}}.$

Since $x^{1},x^{2}\in\mathcal{L}(f(x^{0}))$ assumption (56) yields

$\displaystyle\sup_{f^{*}\in f(\mathcal{L^{*}})}\inf_{z\in
f^{-1}(f^{*})}\frac{2\left(\lVert x^{1}-z\rVert^{2}+\lVert
x^{2}-z\rVert^{2}\right)}{s(k+1)^{2}}\leq\frac{4R}{s(k+1)^{2}},$

with

$\displaystyle R=\max_{j=1,2}\left\\{\sup_{f^{*}\in
f(\mathcal{L^{*}})}\inf_{z\in f^{-1}(f^{*})}\lVert
x^{j}-z\rVert^{2}\right\\}.$ (57)

It remains to show that

$\displaystyle\sup_{z\in\mathcal{H}}\sigma_{k}(z)=\sup_{f^{*}\in
f(\mathcal{L^{*}})}\inf_{z\in f^{-1}(f^{*})}\sigma_{k}(z).$

Writing out the definition of $\sigma_{k}(z)$, we get

$\displaystyle\sup_{f^{*}\in f(\mathcal{L^{*}})}\inf_{z\in
f^{-1}(f^{*})}\sigma_{k}(z)=\sup_{f^{*}\in f(\mathcal{L^{*}})}\inf_{z\in
f^{-1}(f^{*})}\min_{i=1,\dots,m}\left(f_{i}(x^{k})-f_{i}(z)\right)$
$\displaystyle=$ $\displaystyle\sup_{f^{*}\in
f(\mathcal{L^{*}})}\min_{i=1,\dots,m}\left(f_{i}(x^{k})-f_{i}^{*}\right)=\sup_{z\in\mathcal{L}^{*}}\min_{i=1,\dots,m}\left(f_{i}(x^{k})-f_{i}(z)\right)$
$\displaystyle=$
$\displaystyle\sup_{z\in\mathcal{H}}\min_{i=1,\dots,m}\left(f_{i}(x^{k})-f_{i}(z)\right).$

∎

The function
$u_{0}(x)=\sup_{z\in\mathcal{H}}\min_{i=1,\dots,m}f_{i}(x)-f_{i}(z)$ attains
the value zero if and only if $x$ is weakly Pareto optimal. Theorem 6.7 shows
that $u_{0}(x^{k})=\mathcal{O}(k^{-2})$.

###  6.5 Relation to Tanabe’s Accelerated Multiobjective Gradient Method

In the recent preprint [9], Tanabe, Fukuda and Yamashita define an accelerated
proximal gradient method for MOPs with objective functions that have a
separable structure of the form $f_{i}(x)=g_{i}(x)+h_{i}(x)$, where
$g_{i}:\mathbb{R}^{n}\to\mathbb{R}$ is convex, continuously differentiable
with $L$-Lipschitz continuous gradient and $h_{i}:\mathbb{R}^{n}\to\mathbb{R}$
is convex, lower semicontinuous and proper for all $i=1,\dots,m$. Since we
only treat the case of smooth objective functions $f_{i}$, we set from here on
$h_{i}\equiv 0$. Tanabe et al. discovered their method using techniques
different from the ones used throughout this paper, using the concept of merit
functions. We will not recite their method here but refer the reader to [9].
To understand the similarity between their method and Algorithm 1, we
investigate the quadratic optimization problems that have to be solved in each
iteration of the methods, respectively. In the method from [9], the step
direction is computed by solving a quadratic optimization problem with the
following objective function $\Psi:\mathbb{R}^{m}\to\mathbb{R}$,

$\displaystyle\Psi(\theta)\coloneqq\frac{s}{2}\lVert\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k})\rVert^{2}+\sum_{i=1}^{m}\theta_{i}\left(f_{i}(x^{k})-f_{i}(y^{k})\right).$

Using the first order approximation
$f_{i}(y^{k})-f_{i}(x^{k})\approx\langle\nabla
f_{i}(y^{k}),y^{k}-x^{k}\rangle$, we get

$\displaystyle\Psi(\theta)\approx\frac{s}{2}\lVert\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k})\rVert^{2}+\langle\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k}),x^{k}-y^{k}\rangle.$

Minimizing $\Psi(\theta)$ is equivalent to minimizing the function
$\Phi:\mathbb{R}^{m}\to\mathbb{R}$,

$\displaystyle\Phi(\theta)\coloneqq$
$\displaystyle\frac{s^{2}}{2}\lVert\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k})\rVert^{2}+\langle s\sum_{i=1}^{m}\theta_{i}\nabla
f_{i}(y^{k}),x^{k}-y^{k}\rangle+\frac{1}{2}\lVert x^{k}-y^{k}\rVert^{2}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\lVert
s\sum_{i=1}^{m}\theta_{i}\nabla f_{i}(y^{k})+(x^{k}-y^{k})\rVert^{2}.$

Using $x^{k}-y^{k}=-\frac{k-1}{k+2}(x^{k}-x^{k-1})$ we note that
$\Phi(\theta)$ is in fact the objective function of the quadratic optimization
problem (47). After this observation, it is not surprising that the method by
Tanabe et al. shows covergence behavior similar to Algorithm 1.

##  7 Improving the Numerical Efficiency

First order methods for multiobjective optimization that are based on the
steepest descent method by Fliege and Svaiter [12] require the solution of a
quadratic subproblem in each iteration. Computing the solutions of these
problems is computational demanding. In the following subsection, we present a
possible approach to overcome this problem.

###  7.1 A Multiobjective Gradient Method without Quadratic Subproblems

In this subsection, we define a method based on Algorithm 1 which does not
require the solution of a quadratic subproblem in each iteration. In
Subsection 6.2, we derived Algorithm 1 from the scheme

$\displaystyle\frac{3}{k}(x^{k+1}-x^{k})+\operatorname*{proj}_{sC(y^{k})+(x^{k+1}-2x^{k}+x^{k-1})}(0)=0,$
(58)

which can be interpreted as a discretization of the differential equation

$\displaystyle\frac{3}{t}\dot{x}(t)+\operatorname*{proj}_{C(x(t))+\ddot{x}(t)}(0)=0.$

If, instead, we use the discretization

$\displaystyle\frac{3}{k}(x^{k}-x^{k-1})+\operatorname*{proj}_{sC(y^{k})+(x^{k+1}-2x^{k}+x^{k-1})}(0)=0,$
(59)

we obtain a different method. Lemma A.1 gives a formula to compute $x^{k+1}$

$\displaystyle x^{k+1}$
$\displaystyle=-\frac{3}{k}(x^{k}-x^{k-1})-s\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k})+2x^{k}-x^{k-1},$ (60)
$\displaystyle=x^{k}+\frac{k-3}{k}(x^{k}-x^{k-1})-s\sum_{i=1}^{m}\theta_{i}^{k}\nabla
f_{i}(x^{k}),$ (61)

where $\theta^{k}\in\mathbb{R}^{m}$ is a solution to the problem

$\displaystyle\begin{split}\min-\sum_{i=1}^{m}\theta_{i}\langle\nabla
f_{i}(x^{k}),x^{k}-x^{k-1}\rangle\,\,\operatorname{s.t.}\,\,\theta\geq
0\,\,\text{and}\,\,\sum_{i=1}^{m}\theta_{i}=1.\end{split}$ (62)

This can be solved efficiently by computing $m$ inner products. After changing
$\frac{k-3}{k}$ into $\frac{k-1}{k+2}$ in (61), we define Algorithm 2.

Algorithm 2 Accelerated multiobjective gradient method without quadratic
subproblems

1:Choose $x^{0}=x^{1}\in\mathcal{H}$, $s>0$ and set k = 1.

2:Set $y^{k}=x^{k}+\frac{k-1}{k+2}(x^{k}-x^{k-1})$.

3:Compute $j=\operatorname*{arg\,max}_{i=1,\dots,m}\langle\nabla
f_{i}(y^{k}),x^{k}-x^{k-1}\rangle$.

4:Set $x^{k+1}=y^{k}-s\nabla f_{j}(y^{k})$

5:if stopping condition is true then

6: Stop.

7:else

8: Update $k\leftarrow k+1$ and go to step 1.

9:end if

There is no proof of convergence for the method defined by Algorithm 2 but we
discuss its numerical behavior in Section 8. Also, convergence can be
guaranteed by switching from the significantly faster Algorithm 2 to Algorithm
1 as soon as some heuristic criterion is met.

###  7.2 Backtracking for unknown Lipschitz Constants

In all presented algorithms, we can include backtracking if the Lipschitz
constants of the gradients $\nabla f_{i}$ of the objective functions are
unknown. We can do this as stated in [33, 9]. To include backtracking, we
choose an initial step size $s_{0}>0$ and a parameter $\sigma\in(0,1)$. In all
discussed algorithms there is a step $x^{k+1}=w^{k}-sd^{k}$, with
$d^{k}\in\mathcal{H}$ and $w^{k}=x^{k}$ or $w^{k}=y^{k}$. One can replace this
step with $x^{k+1}=w^{k}-s_{k}d^{k}$, with a step size $s_{k}$ that is
determined using backtracking. We choose in every step
$s_{k}=\sigma^{l_{k}}s_{k-1}$ where $l_{k}\geq 0$ is the smallest nonnegative
integer satisfying for all $i=1,\dots,m$

$\displaystyle f_{i}(w^{k}-\sigma^{l_{k}}s_{k-1}d^{k})\leq
f_{i}(w^{k})-\sigma^{l_{k}}s_{k-1}\langle\nabla
f_{i}(w^{k}),d^{k}\rangle+\frac{\sigma^{l_{k}}s_{k-1}}{2}\lVert
d^{k}\rVert^{2}.$

The sequence $(s_{k})_{k\geq 0}$ is monotonically decreasing by definition.
Under the condition that the objective functions posses $L$-Lipschitz
continuous gradients, it is guaranteed that the sequence $(s_{k})_{k\geq 0}$
is constant from same $k$ on. This is true since $s_{k}$ can only decrease as
long as $s_{k}>\frac{1}{L}$. Therefore, $s_{k}$ can only decrease finitely
many times until it reaches a point where $s_{k}\leq\frac{1}{L}$. Using this
observation, we can include backtracking in Algorithm 1 and still use the
proofs of Theorem 6.6 and Theorem 6.7 to show that the same convergence
results can be achieved.

##  8 Numerical Examples

In this section, we present the typical behavior of our algorithms on two test
problems. We compare Algorithms 1 and 2 with the steepest descent method by
Fliege and Svaiter with constant step sizes [12]. Throughout this section we
denote the steepest descent method by SD, Algorithm 1 by AccG (accelerated
gradient method) and Algorithm 2 by AccG w\o Q (accelerated gradient method
without quadratic subproblems). We implemented all codes using Matlab R2021b
and executed the algorithms on a machine with a 2.80 GHz Intel Core i7
processor and 48 GB memory. We solved the quadratic subproblems for SD and
AccG using the built-in Matlab function quadprog.

###  8.1 Example 1: A Convex MOP with three Objective Functions

(a)

(b)

Figure 1: Coordinates $(x_{1},x_{2},x_{3})$ of the sequences $(x^{k})_{k\geq
0}$ for SD, AccG and AccG w\o Q. Line plot for $1000$ iterations with a filled
circle every $50$ iterations to compare the velocities.

(c) Function values $f_{1}(x^{k})$

(d) Function values $f_{2}(x^{k})$

(e) Function values $f_{3}(x^{k})$

Figure 2: Function values $(f_{i}(x^{k}))_{k\geq 0}$ of the iterates for the
objective functions $i=1,2,3$ for the different algorithms.

(a) $k_{\max}=50$

(b) $k_{\max}=250$

(c) $k_{\max}=1000$

(d) $k_{\max}=1000$

Figure 3: Function values of the objective functions in the image space for
the different algorithms and a different maximum number of iterations
$k_{\max}=50,250,1000$.

In our first example, we choose a problem with input dimension $n=20$ and
three objective functions ($m=3$). We define the objective functions using the
following parameters. For $p=50$ and $i=1,2,3$ we generate matrices
$A^{i}=\left(a_{1}^{i},\dots,a_{p}^{i}\right)^{T}\in\mathbb{R}^{p\times n}$
with $a_{j}^{i}\in\mathbb{R}^{n}$ for $j=1,\dots,p$ and vectors
$b^{i}\in\mathbb{R}^{p}$. Then, for $i=1,2,3$, we define the objective
functions

$\displaystyle f_{i}:\mathbb{R}^{n}\to\mathbb{R},\quad
x\mapsto\ln\left(\sum_{j=1}^{p}\exp\left((a_{j}^{i})^{T}x\right)\right).$

For the first experiment we randomly generate Matrices
$A^{i}\in\mathbb{R}^{p\times n}$ and vectors $b_{i}\in\mathbb{R}^{p}$ with
entries uniformly sampled in $[-1,1]$ for $i=1,2,3$. The starting vector
$x_{0}$ is uniformly randomly drawn from $[-15,15]^{n}$. We use the step size
$s={5}\mathrm{e}{-2}$ and execute maximally $k_{\max}=1000$ iterations. Figure
2 contains plots of the sequences $(x^{k})_{k\geq 0}$ for the different
algorithms. In Figure 1(a), one sees that the sequences generated with AccG
and AccG w\o Q advance much faster in the beginning, while the velocity for
the sequence generated with SD remains constant. The sequences generated by
AccG and AccG w\o Q give very similar trajectories in the beginning. This
result is intuitive given that the schemes in the algorithms are derived from
different discretizations of the same differential equation. However, this
result is still surprising keeping in mind that in Algorithm 2 we do not solve
a quadratic subproblem in each iteration. Only in Figure 1(b), we see that the
sequences differ more substantially in the long run. It is also noteworthy
that the sequence generated by AccG is smoother compared to the trajectory
generated by AccG w\o Q. This is due to the fact that in AccG w\o Q we choose
one of the gradients of the objective functions for the gradient component of
the step direction while in AccG we choose an element of the convex hull of
the gradients. AccG and AccG w\o Q are superior to SD in terms of convergence
of the function values for all objective functions, as shown in Figure 2. AccG
and AccG w\o Q experience fast convergence within the first $200$ iterations.
Comparing the different objective functions in Figures 1(c), 1(d) and 1(e), we
see that AccG and AccG w\o Q yield outputs with similar function values for
all objective functions.

In a second experiment, we execute all algorithms for $50$ starting values
uniformly sampled in $[-5,5]^{n}$ with step size $s={5}\mathrm{e}{-2}$. We use
the stopping criterion $\lVert
f(x^{k})-f(x^{k-1})\rVert_{\infty}<{1}\mathrm{e}{-4}$ to stop the algorithms
if the function values do not change significantly. In Figure 3(a), we perform
up to $k_{\max}=50$, in Figure 3(b) up to $k_{\max}=250$ and in Figures 3(c)
and 3(d) up to $k_{\max}=1000$ iterations. Similar to the results observed in
Figures 2 and 2, AccG and AccG w\o Q advance much faster in the beginning
compared to SD. Comparing Figures 1(d) and 1(e), we see that after $250$
iterations the function values for the accelerated methods are converging or
the stopping conditions were met. The different behavior of the accelerated
methods can be observed in Figure 3(d). While the solutions of AccG are
farther spread, it looks like the solutions of AccG w\o Q are drawn towards
the center of the Pareto front. Altogether, the accelerated methods perform
better for this problem in terms of convergence speed of the function values.
In Table 1 the total number of iterations and computation times for the
experiments are listed. The accelerated methods require fewer iterations.
Compared to SD, AccG requires only approximately 25 % and AccG w\o Q only
approximately 50 % iterations. For the computation times the results are
different. SD and AccG behave similar, with AccG requiring approximately 25 %
of the computation time that is required for SD. However, AccG w\o Q needs
less the 2 % of the time which is consumed by AccG. This improvement stems
from the quadratic optimization problems that not required in AccG w\o Q.

| SD | AccG | AccG w\o Q
---|---|---|---
total iterations | 49924 | 12230 | 25906
total time | $436.54\,\text{s}$ | $100.94\,\text{s}$ | $1.61\,\text{s}$
Table 1: Total iterations and computation times for algorithm executions using
parameters $s={5}\mathrm{e}{-2}$, $k_{\max}=1000$ and stopping condition
$\lVert f(x^{k})-f(x^{k-1})\rVert_{\infty}<{1}\mathrm{e}{-4}$ for $50$ start
values uniformly sampled in $[-15,15]^{n}$.

###  8.2 Example 2: A Nonconvex MOP with two Objective Functions

For our second test problem, we choose an example from [34] with input
dimension $n=2$ and the two objective functions

$\displaystyle
f_{1}(x)=\frac{1}{2}(\sqrt{1+(x_{1}+x_{2})^{2}}+\sqrt{1+(x_{1}-x_{2})^{2}}+x_{1}-x_{2})+\lambda\exp{(-(x_{1}-x_{2})^{2})},$
$\displaystyle
f_{2}(x)=\frac{1}{2}(\sqrt{1+(x_{1}+x_{2})^{2}}+\sqrt{1+(x_{1}-x_{2})^{2}}-x_{1}+x_{2})+\lambda\exp{(-(x_{1}-x_{2})^{2})},$

with $\lambda=0.6$. For the multiobjective optimization problem (8) with these
objective functions it can easily be verified that the Pareto set is

$\displaystyle P=\left\\{x\in\mathbb{R}^{2}\,\,:\,\,x_{1}+x_{2}=0\right\\}.$

In the first experiment, we execute Algorithms SD, AccG and AccG w\o Q with
the starting vector $x^{0}=(1,2)^{T}$ and perform $k_{\max}=1000$ iterations.
The step size is set to $s={5}\mathrm{e}{-3}$. The sequences and function
values of the objective functions are shown in Figure 4. Similarly to the
first experiment in Figure 4(a), the sequences $(x^{k})_{k\geq 0}$ of the
accelerated methods advance faster in the beginning. While Algorithms SD and
AccG converge to the same element in the Pareto set the algorithm AccG w\o Q
produces a trajectory that deviates from the trajectories of SD and AccG and
moves to a different part of the Pareto set. The values of the objective
functions in Figures 4(b) and 4(c) indicate a similar behavior. For the
accelerated methods we have faster decrease in the beginning and we note that
the function values for SD and AccG converge to similar values. In Figure 5 we
use $100$ random starting points uniformly sampled in $[-2,2]^{2}$. For the
experiments we use different maximal numbers of iterations $k_{\max}$. In
addition we stop the algorithm if $\lVert
f(x^{k})-f(x^{k-1})\rVert_{\infty}<{1}\mathrm{e}{-4}$. Comparing Figures 5(a),
5(b) and 5(c), we note that the objective function values of the accelerated
methods decrease much faster in the beginning. For $k_{\max}=100$ Algorithms
AccG and AccG w\o Q yield solutions that are distributed along the Pareto
front. In Table 2 we list the total number of iterations and total computation
times for executions with up to $k_{\max}=1000$ iterations with stopping
condition $\lVert f(x^{k})-f(x^{k-1})\rVert_{\infty}<{1}\mathrm{e}{-4}$.
Compared to SD, AccG needs only approximately 15 % and AccG w\o Q only
approximately 51 % iterations.

(a) Sequences $(x^{k})_{k\geq 0}$

(b) Function values $f_{1}(x^{k})$

(c) Function values $f_{2}(x^{k})$

Figure 4: Sequences $(x^{k})_{k\geq 0}$ and function values
$(f_{i}(x^{k}))_{k\geq 0}$ of iterates for $i=1,2$. For the sequences we use a
line plot for 1000 iterations with a filled circle every $50$ iterations to
compare velocities.

(a) $k_{\max}=50$

(b) $k_{\max}=100$

(c) $k_{\max}=500$

Figure 5: Values of the objective functions in the image space for different maximum numbers of iterations $k_{\max}=50,250,1000$ for the different algorithms SD, AccG and AccG w\o Q | SD | AccG | AccG w\o Q
---|---|---|---
total iterations | 45543 | 6632 | 23034
total time | $431.75\,\text{s}$ | $62.90\,\text{s}$ | $0.31\,\text{s}$
Table 2: Total iterations and computation times for algorithm executions using
parameters $h={5}\mathrm{e}{-3}$, $k_{\max}=1000$ and stopping condition
$\lVert f(x^{k})-f(x^{k-1})\rVert_{\infty}<{1}\mathrm{e}{-4}$ for $100$ start
values uniformly sampled in $[-2,2]^{n}$.

In another experiment we compare how the choice of the step size $s$ affects
the solutions of AccG w\o Q. We use the step sizes
$s={5}\mathrm{e}{-3},{1}\mathrm{e}{-2},{5}\mathrm{e}{-2}$. For all executions
we perform $k_{\max}=1000$ iterations, with the stopping criterion $\lVert
f(x^{k})-f(x^{k-1})\rVert_{\infty}<{1}\mathrm{e}{-4}$. Comparing Figures 6(a),
6(b) and 6(c) we see that for the smallest step size $s={5}\mathrm{e}{-3}$
solutions are distributed on the whole Pareto front. For the biggest step size
$s={5}\mathrm{e}{-2}$ Algorithm AccG w\o Q yields solutions that cluster at
two points of the Pareto front, which is not desirable in general. The two
points where the solutions cluster correspond to the knee points in the Pareto
front. This is not surprising since these points correspond to solutions where
the individual gradients are similar in magnitude, which is why the solution
is zig-zagging back and forth between the objectives in these locations.
However, this disadvantage is compensated by the fact that we do not need to
solve a quadratic subproblem in every step. We need potentially more
iterations when choosing smaller step sizes but every iteration is
computationally cheaper in comparison to SD and AccG. Moreover, a small
adaptation to step 2 of Algorithm 2, where we include a weighting parameter in
the $\max$ problem, might allow us to diversify solutions, similar to the
weighted sum method.

(a) $s={5}\mathrm{e}{-3}$

(b) $s={1}\mathrm{e}{-2}$

(c) $s={5}\mathrm{e}{-2}$

Figure 6: Values of the objective functions in the image space for different
step sizes $s={5}\mathrm{e}{-3},{1}\mathrm{e}{-2},{5}\mathrm{e}{-2}$ for the
algorithm AccG w\o Q.

##  9 Conclusion and Open Questions

We present the novel inertial gradient-like dynamical system (IMOG’) for
Pareto optimization. We show that trajectories of this system converge weakly
to Pareto critical points of (8). Based on this, we define a novel inertial
gradient method for multiobjective optimization and show weak convergence to
Pareto critical points. We derive an accelerated gradient method from the
informally introduced internal gradient-like system (MAVD) which incorporates
asymptotically vanishing damping. Using the concept of merit functions, we
show that our method possesses an improved convergence rate. Using a different
discretization of the system (MAVD), we define a further accelerated gradient
method which does not require the solution to a quadratic optimization problem
in every step. A comparison on selected test problems shows that the
accelerated methods are in fact superior to the plain multiobjective steepest
descent method.

There are a lot of open questions arising from the presented work. The
gradient system (IMOG’) can be analyzed for different problem classes. Also,
we have not formally discussed the system (MAVD). It would be interesting to
see if we can derive solutions from this system for $\alpha>3$ and obtain
stronger convergence results similar to the single objective setting. In
addition, we can adapt our gradient systems and algorithms to treat problems
with a separable smooth and nonsmooth structure using proximal methods. It
would also be interesting to investigate the behavior of our algorithms for
high-dimensional and nonconvex problems. In addition, one could apply the
presented algorithms in the area of machine learning, e.g., for multitask
learning problems [35].

##  Appendix A Two Lemmas on Convex Projections

######  Lemma A.1.

Let $\mathcal{H}$ be a real Hilbert space, $C\subset\mathcal{H}$ a convex and
compact set and $\eta\in\mathcal{H}$ a fixed vector. Then, $\xi\in\mathcal{H}$
is a solution to the problem

$\displaystyle\text{Find $\xi\in\mathcal{H}$ such that :
}\,\eta=\operatorname*{proj}_{C+\xi}(0),$ (63)

if and only if it has the form $\xi=\eta-\mu$, where $\mu$ is a solution to
the constrained optimization problem $\min_{\mu\in
C}\,\langle\mu,\eta\rangle$.

###### Proof.

First, we show that an element of the form $\xi=\eta-\mu$, with $\mu$ a
solution to $\min_{\mu\in C}\,\langle\mu,\eta\rangle$ is a solution to problem
(63). The set of minimizers of the problem $\min_{\mu\in
C}\,\langle\mu,\eta\rangle$ is nonempty, since $C$ is compact. Let
$\mu\in\operatorname*{arg\,min}_{\mu\in C}\langle\mu,\eta\rangle$. Since $C$
is convex, the first order optimality condition for this problem gives that
for all $x\in C$ it holds that $\langle x-\mu,\eta\rangle\geq 0$ and hence

$\displaystyle\langle x+\xi-(\mu+\xi),\eta\rangle\geq 0.$

Since we have chosen $\xi=\eta-\mu$ the equation above reads as

$\displaystyle\langle x+\xi-\eta,\eta\rangle\geq 0,$

which is equivalent to $\eta=\operatorname*{proj}_{C+\xi}(0)$. The other
direction works analogously. If the vector $\xi$ is a solution to problem (63)
this guarantees that $\mu=\xi-\eta$ satisfies the first order optimality
condition for problem $\min_{\mu\in C}\,\langle\mu,\eta\rangle$. Since problem
$\min_{\mu\in C}\,\langle\mu,\eta\rangle$ is convex and defined over a convex
set, this is equivalent to $\mu$ being an optimal solution to $\min_{\mu\in
C}\,\langle\mu,\eta\rangle$. ∎

######  Lemma A.2.

Let $\mathcal{H}$ be a real Hilbert space, $C\subset\mathcal{H}$ a convex and
closed set and $a>0,\nu\in\mathcal{H}$ fixed, with $a\neq-1$ . Then, the
problem

$\displaystyle\text{Find $\xi\in\mathcal{H}$ such that :
}\,-a(\xi+\nu)=\operatorname*{proj}_{C+\xi}(0),$ (64)

has the unique solution
$\xi=-\left(\frac{1}{1+a}\operatorname*{proj}_{C}\left(\nu\right)+\frac{a}{1+a}\nu\right)$.

###### Proof.

First, we show that
$\xi=-\left(\frac{1}{1+a}\operatorname*{proj}_{C}\left(\nu\right)+\frac{a}{1+a}\nu\right)$
is a solution to (64). It is easy to check that $-a(\xi+\nu)\in C+\xi$. Define
$p\coloneqq\operatorname*{proj}_{C}\left(\nu\right)$. For all $x\in C$ it
holds that $\langle x-p,p-\nu\rangle\geq 0$ and hence for all $x\in C$ we get

$\displaystyle\langle x+\xi+a(\xi+\nu),a(\xi+\nu)\rangle\leq 0,$

which is equivalent to

$\displaystyle-a(\xi+\nu)=\operatorname*{proj}_{C+\xi}(0).$

The uniqueness follows the same way. Assume we have a solution $\tilde{\xi}$
to (64). By the same computations as above it holds that for all $x\in C$

$\displaystyle\langle x+(1+a)\tilde{\xi}+a\nu,\tilde{\xi}+\nu\rangle\leq 0.$

This is equivalent to

$\displaystyle-((1+a)\tilde{\xi}+a\nu)=\operatorname*{proj}_{C}(\nu),$

from which follows that $\xi=\tilde{\xi}$ is the unique solution. ∎

## References

  * [1] Weijie Su, Stephen Boyd, and Emmanuel Candes.  A differential equation for modeling nesterov’s accelerated gradient method: theory and insights.  Advances in neural information processing systems, 27, 2014.
  * [2] Hedy Attouch and Guillaume Garrigos.  Multiobjective optimization: an inertial dynamical approach to pareto optima.  arXiv preprint arXiv:1506.02823, 2015.
  * [3] Boris T Polyak.  Some methods of speeding up the convergence of iteration methods.  Ussr computational mathematics and mathematical physics, 4(5):1–17, 1964.
  * [4] Felipe Alvarez.  On the minimizing property of a second order dissipative system in hilbert spaces.  SIAM Journal on Control and Optimization, 38(4):1102–1119, 2000\.
  * [5] Hedy Attouch, Xavier Goudou, and Patrick Redont.  The heavy ball with friction method, i. the continuous dynamical system: global exploration of the local minima of a real-valued function by asymptotic analysis of a dissipative dynamical system.  Communications in Contemporary Mathematics, 2(01):1–34, 2000.
  * [6] Hedy Attouch and Juan Peypouquet.  The rate of convergence of nesterov’s accelerated forward-backward method is actually $o(k^{-2})$.  SIAM Journal on Optimization, 26, 10 2015.
  * [7] Hedy Attouch, Zaki Chbani, Juan Peypouquet, and Patrick Redont.  Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity.  Mathematical Programming, 168(1):123–175, 2018.
  * [8] Hedy Attouch, Zaki Chbani, and Hassan Riahi.  Rate of convergence of the nesterov accelerated gradient method in the subcritical case $\alpha\leq 3$.  ESAIM: Control, Optimisation and Calculus of Variations, 25:2, 2019\.
  * [9] Hiroki Tanabe, Ellen H Fukuda, and Nobuo Yamashita.  An accelerated proximal gradient method for multiobjective optimization.  arXiv preprint arXiv:2202.10994, 2022.
  * [10] Kaisa Miettinen.  Nonlinear Multiobjective Optimization.  Springer, New York, 1998.
  * [11] Mustapha El Moudden and Abdelkrim El Mouatasim.  Accelerated diagonal steepest descent method for unconstrained multiobjective optimization.  Journal of Optimization Theory and Applications, 188(1):220–242, 2021.
  * [12] Jörg Fliege and Benar Fux Svaiter.  Steepest descent methods for multicriteria optimization.  Mathematical methods of operations research, 51(3):479–494, 2000\.
  * [13] Yurii E Nesterov.  A method for solving the convex programming problem with convergence rate o (1/k^ 2).  In Dokl. akad. nauk Sssr, volume 269, pages 543–547, 1983.
  * [14] Steve Smale.  Global analysis and economics i: Pareto optimum and a generalization of morse theory.  In Dynamical systems, pages 531–544. Elsevier, 1973.
  * [15] Claude Henry.  An existence theorem for a class of differential equations with multivalued right-hand side.  Journal of Mathematical Analysis and Applications, 41(1):179–186, 1973.
  * [16] Bernard Cornet.  Monotone planning procedures and accessibility of pareto optima.  New trends in dynamic system theory and economics. Academic Press, Cambridge, pages 337–349, 1979.
  * [17] B Cornet.  Contributions to the mathematical theory of dynamical mechanisms of resource allocation.  These de doctorat d’état, Université Paris IX Dauphine, 1981.
  * [18] Bernard Cornet.  Existence of slow solutions for a class of differential inclusions.  Journal of mathematical analysis and applications, 96(1):130–147, 1983.
  * [19] Enrico Miglierina.  Slow solutions of a differential inclusion and vector optimization.  Set-Valued Analysis, 12(3):345–356, 2004.
  * [20] Hedy Attouch and Xavier Goudou.  A continuous gradient-like dynamical approach to pareto-optimization in hilbert spaces.  Set-Valued and Variational Analysis, 22(1):189–219, 2014.
  * [21] Hedy Attouch, Guillaume Garrigos, and Xavier Goudou.  A dynamic gradient approach to pareto optimization with nonsmooth convex objective functions.  Journal of Mathematical Analysis and Applications, 422(1):741–771, 2015.
  * [22] Hiroki Tanabe, Ellen H. Fukuda, and Nobuo Yamashita.  New merit functions for multiobjective optimization and their properties, 2020.
  * [23] Guang-Ya Chen, Chuen-Jin Goh, and Xiao Qi Yang.  On gap functions for vector variational inequalities.  In Vector variational inequalities and vector equilibria, pages 55–72. Springer, 2000.
  * [24] X Q Yang and J C Yao.  Gap functions and existence of solutions to set-valued vector variational inequalities.  Journal of Optimization Theory and Applications, 115(2):407–417, 2002.
  * [25] C G Liu, K F Ng, and W H Yang.  Merit functions in vector optimization.  Mathematical programming, 119(2):215–237, 2009.
  * [26] J-P Aubin and Arrigo Cellina.  Differential inclusions: set-valued maps and viability theory, volume 264.  Springer Science & Business Media, 2012.
  * [27] Haim Brezis.  Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert.  Elsevier, 1973.
  * [28] Zdzisław Opial.  Weak convergence of the sequence of successive approximations for nonexpansive mappings.  Bulletin of the American Mathematical Society, 73(4):591–597, 1967\.
  * [29] Felipe Alvarez and Hedy Attouch.  An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping.  Set-Valued Analysis, 9:3–11, 01 2001.
  * [30] Heda Attouch and Jalal Fadili.  From the ravine method to the nesterov method and vice versa: a dynamical system perspective.  arXiv preprint arXiv:2201.11643, 2022.
  * [31] Hiroki Tanabe, Ellen H Fukuda, and Nobuo Yamashita.  Convergence rates analysis of a multiobjective proximal gradient method.  Optimization Letters, pages 1–18, 2022.
  * [32] David G Luenberger.  Optimization by vector space methods.  John Wiley & Sons, 1997.
  * [33] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM journal on imaging sciences, 2(1):183–202, 2009.
  * [34] Katrin Witting.  Numerical algorithms for the treatment of parametric multiobjective optimization problems and applications.  PhD thesis, Paderborn, Universität Paderborn, Diss., 2012, 2012.
  * [35] O. Sener and V. Koltun.  Multi-task learning as multi-objective optimization.  In Advances in Neural Information Processing Systems, volume 31, 2018.
  * [36] Felipe Alvarez, Hedy Attouch, Jérôme Bolte, and Patrick Redont.  A second-order gradient-like dissipative dynamical system with hessian-driven damping.: Application to optimization and mechanics.  Journal de mathématiques pures et appliquées, 81(8):747–779, 2002.
  * [37] Haim Brezis.  Monotonicity methods in hilbert spaces and some applications to nonlinear partial differential equations.  In Contributions to nonlinear functional analysis, pages 101–156. Elsevier, 1971.
  * [38] Ronald E Bruck Jr.  Asymptotic convergence of nonlinear contraction semigroups in hilbert space.  Journal of Functional Analysis, 18(1):15–26, 1975.
  * [39] Jörg Fliege, A Ismael F Vaz, and Luís Nunes Vicente.  Complexity of gradient descent for multiobjective optimization.  Optimization Methods and Software, 34(5):949–959, 2019.
  * [40] Hedy Attouch and Alexandre Cabot.  Convergence rates of inertial forward-backward algorithms.  SIAM Journal on Optimization, 28(1):849–874, 2018.
  * [41] Bennet Gebken and Sebastian Peitz.  An efficient descent method for locally lipschitz multiobjective optimization problems.  Journal of Optimization Theory and Applications, 188(3):696–723, 2021.
  * [42] Joerg Fliege, L M Grana Drummond, and Benar Fux Svaiter.  Newton’s method for multiobjective optimization.  SIAM Journal on Optimization, 20(2):602–626, 2009.
  * [43] Michael Dellnitz, Oliver Schütze, and Thorsten Hestermeyer.  Covering pareto sets by multilevel subdivision techniques.  Journal of optimization theory and applications, 124(1):113–136, 2005.
  * [44] Hedy Attouch, Radu Ioan Boţ, and Ernö Robert Csetnek.  Fast optimization via inertial dynamics with closed-loop damping.  Journal of the European Mathematical Society, 2022.
  * [45] Sebastian Peitz, Sina Ober-Blöbaum, and Michael Dellnitz.  Multiobjective optimal control methods for the navier-stokes equations using reduced order modeling.  Acta Applicandae Mathematicae, 161(1):171–199, 2019.
  * [46] Sina Ober-Blöbaum and Sebastian Peitz.  Explicit multiobjective model predictive control for nonlinear systems with symmetries.  International Journal of Robust and Nonlinear Control, 31(2):380–403, 2021.