# An effective subgradient algorithm via Mifflin’s line search for nonsmooth
nonconvex multiobjective optimization

Morteza Maleknia  School of Mathematics, Statistics and Computer Science,
College of Science, University of Tehran, Tehran, Iran  Majid Soleimani-
damaneh  Department of Mathematical Sciences, Isfahan University of
Technology, Isfahan, Iran

## Abstract

We propose a descent subgradient algorithm for unconstrained nonsmooth
nonconvex multiobjective optimization problems. To find a descent direction,
we present an iterative process that efficiently approximates the (Goldstein)
$\varepsilon$-subdifferential of each objective function. To this end, we
develop a new variant of Mifflin’s line search in which the subgradients are
arbitrary and its finite convergence is proved under a semismooth assumption.
In order to reduce the number of subgradient evaluations, we employ a
backtracking line search that identifies the objectives requiring an
improvement in the current approximation of the $\varepsilon$-subdifferential.
Meanwhile, for the remaining objectives, new subgradients are not computed.
Unlike bundle-type methods, the proposed approach can handle nonconvexity
without the need for algorithmic or parametric adjustments. Moreover, the
quadratic subproblems have a simple structure, and hence the method is easy to
implement. We analyze the global convergence of the proposed method and prove
that any accumulation point of the generated sequence satisfies a necessary
Pareto optimality condition. Furthermore, our convergence analysis addresses a
theoretical challenge in a recently developed subgradient method. By
conducting numerical experiments, we observe the practical capability of the
proposed method and evaluate its efficiency when applied to a diverse range of
nonsmooth test problems.

##  1 Introduction

In many practical situations, one needs to make a decision under multiple
conflicting criteria. For instance, consider a refinery that intends to
maximize petrol production while minimizing toxic emissions at the same time.
This problem is a decision-making problem with two conflicting objectives.
Such a problem is called a _M ultiobjective Optimization Problem_ (MOP). In
contrast with single objective problems, due to the conflicts among
objectives, there is no solution to an MOP optimizing all objectives
simultaneously. Therefore, the optimality notion is replaced with _Pareto
optimality_ , which represents an optimal compromise among all objectives. In
practical terms, it is of crucial importance to approximate the set of all
Pareto optimal points, namely the _Pareto front_ , which is a challenging
issue. This problem becomes even more complicated when objectives are not
smooth and convex.

In this paper, we consider the following unconstrained multiobjective
minimization problem:

$\min\,\,\left(f_{1}(\mathbf{x}),f_{2}(\mathbf{x}),\ldots,f_{p}(\mathbf{x})\right)\quad\text{s.t.}\quad\mathbf{x}\in\mathbb{R}^{n},$
(1)

where the objective functions $f_{i}:\mathbb{R}^{n}\to\mathbb{R}$,
$i=1,\ldots,p$, are locally Lipschitz, but not necessarily smooth or convex.
One of the main approaches to deal with MOP (1) is to employ the scalarization
methods [1, 2, 3]. This approach turns an MOP to a nonsmooth single-objective
optimization problem. Then, a global optimal solution of the resulted single-
objective problem provides a (weak) Pareto optimal point for the MOP.
Regarding a nonsmooth single-objective optimization problem, various types of
methods are developing rapidly. For instance, subgradient methods [4, 5],
bundle methods [6, 7, 8, 9, 10], and gradient sampling methods [11, 12, 13,
14] are among the most common numerical approaches to minimize a locally
Lipschitz function. The weighted sum method [15, 16] is widely recognized as
one of the most popular scalarization methods. By adjusting the weighting
vectors, it is possible to obtain the entire Pareto front in the convex case,
at least in theory. However, in situations where the image space is not
$\mathbb{R}^{p}_{\geq}$-convex, which occurs in nonconvex cases, a significant
portion of the Pareto front cannot be achieved using this method (for a bi-
objective example, see Subsection 5.2). In addition, for some MOPs, even in
the convex case, a broad range of weighting vectors lead to an unbounded
single-objective problem (for a bi-objective instance, see [17]). Similar
limitations can also be observed with other aggregation techniques.

As another class of methods for solving MOP (1), one can point to the
iterative descent methods. These approaches are extensions of the existing
descent methods in single-objective optimization to the multiobjective case,
which can be directly applied to the non-scalarized MOP (1). These methods are
descent in the sense that, during each iteration, a descent direction that
simultaneously reduces the value of all the objectives is found and utilized
to update the current point. Due to the theoretical and practical challenges
in developing descent methods for MOPs, the literature along this line of
research is scarce, especially for problems with nonsmooth and nonconvex data.
Some extensions of the subgradient methods that can be applied to MOPs with
convex objectives can be found in [18, 19, 20]. As a more recent subgradient
method, Gebken and Peitz [21] developed a subgradient method that can
efficiently solve nonsmooth nonconvex MOP (1). In [22], Mäkelä et al.
presented an effective generalization of the proximal bundle methods that can
deal with MOP (1). A quasi-Newton method can be found in [23]. Furthermore,
some variants of the proximal bundle methods for solving nonsmooth nonconvex
MOPs which are able to handle inexact information have been proposed in [24,
25].

In this paper, we develop a descent implementable subgradient algorithm to
find Pareto optimal points of MOP (1). The main idea is to develop an
iterative approach to efficiently approximate the (Goldstein)
$\varepsilon$-subdifferential [26] of each objective function. The least-norm
element in the convex hull of the approximated $\varepsilon$-subdifferentials
provides a trial search direction. Such an element is computed by solving a
convex quadratic subproblem. Next, a limited backtracking line search
determines whether we can take a serious step or we need to take a null step
to improve our approximations of $\varepsilon$-subdifferentials. In this way,
for the given positive sequences $\varepsilon_{\nu}\downarrow 0$ and
$\delta_{\nu}\downarrow 0$, we generate a sequence of
$\big{(}\delta_{\nu},\mathcal{S}_{\varepsilon_{\nu}}(\mathbf{x}_{\nu+1})\big{)}$-substationary
points (see Definition 4.1), such that any accumulation point of this sequence
is Clarke substationary for MOP (1).

Regarding the details of the presented algorithm, the following features are
highlighted. The _improvement function_ which is used in Mäkelä et al’s
proximal bundle method [22, 1] poses an underlying issue. Indeed, once a null
step occurs, the improvement function is not able to identify those objectives
whose bundles of subgradient have not been enriched adequately. Therefore,
after a null step, it is inevitable to compute a new subgradient for all the
objectives. The same happens in recently developed proximal bundle methods
[24, 25]. To resolve this difficulty, we propose a simple backtracking line
search that excludes the objectives for which, after a null step, we do not
need to enrich the bundle of information. For the remaining objectives, we
develop a new variant of Mifflin’s line search [6] that finds a new effective
subgradient to improve the bundle of information. In this manner, one can
significantly reduce the number of subgradient evaluations. Moreover, in
Gebken and Peitz [21], the proposed line search procedure to find effective
subgradients lacks the finite convergence property. More precisely, they
merely argue that the asymptotic behavior of the proposed line search provides
a subdifferential set within which the desired subgradient element exists (see
Lemma 3.4 in [21]). In contrast, our proposed line search exactly finds the
desired subgradient element in finite number of iterations.

We study the global convergence of the proposed method and show that any
accumulation point of the generated sequence is Clarke substationary for MOP
(1). Furthermore, in Gebken and Peitz [21], convergence of the method to a
Clarke substationary point through a sequence of
$\big{(}\delta,\mathcal{S}_{\varepsilon}(\mathbf{x})\big{)}$-substationary
points remains as a theoretical challenge (see Remark 4.1 in [21]). This
problem is comprehensively solved in our convergence analysis.

By means of numerical tests, we show the efficiency and ability of the
proposed method in practice. For this purpose, we first analyze the typical
behavior of the method on a bi-objective test problem. Next, by using the
multi-start technique, we assess the efficiency of our method in approximating
the Pareto front. In our third experiment, we seek a sparse solution to an
underdetermined linear system by minimizing a nonsmooth bi-objective model.
Eventually, a set of nonsmooth convex and nonconvex test problems are
considered to provide some comparative results with the multiobjective
proximal bundle method proposed by Mäkelä et al. [27].

The rest of the paper is organized as follows. Some required basic concepts
are provided in Section 1. Our main strategy for finding a descent direction
along with finite convergence property of the proposed line search is given in
Section 3. Global convergence is studied in Section 4. Numerical results are
reported in Section 5, and Section 6 concludes the paper.

##  2 Mathematical preliminaries

Throughout this paper, we use the following notations. The usual inner product
in the Euclidean space $\mathbb{R}^{n}$ is denoted by
$\mathbf{x}^{T}\mathbf{y}$, which induces the Euclidean norm
$\lVert\mathbf{x}\rVert=(\mathbf{x}^{T}\mathbf{x})^{1/2}$. An open ball with
center $\mathbf{x}\in\mathbb{R}^{n}$ and radius $\varepsilon>0$ is denoted by
$\mathcal{B}(\mathbf{x},\varepsilon)$, that is,

$\mathcal{B}(\mathbf{x},\varepsilon):=\\{\mathbf{y}\in\mathbb{R}^{n}\,\,:\,\,\lVert\mathbf{y}-\mathbf{x}\rVert<\varepsilon\\}.$

Moreover, $\mathbb{N}$ is the set of natural numbers,
$\mathbb{N}_{0}:=\mathbb{N}\cup\\{0\\}$, and $\mathbb{R}_{+}:=(0,\infty).$

Suppose $f:\mathbb{R}^{n}\to\mathbb{R}$ is a locally Lipschitz function. Then,
by Rademacher’s theorem [28], $f$ is differentiable almost everywhere on
$\mathbb{R}^{n}$. Let

$\Omega_{f}:=\\{\mathbf{x}\in\mathbb{R}^{n}\,\,:\,\,f\,\,\text{is not
differentiable at}\,\,\mathbf{x}\\}.$

Then, the Clarke subdifferential of $f$ at a point
$\mathbf{x}\in\mathbb{R}^{n}$ is defined as [29]

$\partial
f(\mathbf{x}):=\texttt{conv}\\{\boldsymbol{\xi}\in\mathbb{R}^{n}\,\,:\,\,\exists\,\\{\mathbf{x}_{i}\\}\subset\mathbb{R}^{n}\setminus\Omega_{f}\,\,\,\,\text{s.t.}\,\,\,\,\mathbf{x}_{i}\to\mathbf{x}\,\,\text{and}\,\,\nabla
f(\mathbf{x}_{i})\to\boldsymbol{\xi}\\},$

where conv denotes the convex hull of a set. Furthermore, for any
$\varepsilon>0,$ the (Goldstein) $\varepsilon$-subdifferential of $f$ at a
point $\mathbf{x}\in\mathbb{R}^{n}$ is the set [26]

$\partial_{\varepsilon}f(\mathbf{x}):=\texttt{cl\,conv}\\{\partial
f(\mathbf{y})\,\,:\,\,\mathbf{y}\in\mathcal{B}(\mathbf{x},\varepsilon)\\},$

in which cl conv is the closure of the convex hull. For any $\varepsilon>0$
and $\mathbf{x}\in\mathbb{R}^{n}$, the set
$\partial_{\varepsilon}f(\mathbf{x})$ is a nonempty, convex and compact subset
of $\mathbb{R}^{n}$. If $f$ is differentiable at
$\mathbf{x}\in\mathbb{R}^{n}$, then $\nabla f(\mathbf{x})\in\partial
f(\mathbf{x})$. Furthermore, if $f$ is smooth at
$\mathbf{x}\in\mathbb{R}^{n}$, we have $\partial f(\mathbf{x})=\\{\nabla
f(\mathbf{x})\\}$. Also, the set-valued map
$\partial_{\varepsilon}f:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is
locally bounded and upper semicontinuous [29]. The generalized directional
derivative of function $f$ at point $\mathbf{x}\in\mathbb{R}^{n}$ and
direction $\mathbf{d}\in\mathbb{R}^{n}$ is given by [29]

$f^{\circ}(\mathbf{x};\mathbf{d}):=\limsup_{\begin{subarray}{c}\mathbf{y}\to\mathbf{x}\\\
t\downarrow
0\end{subarray}}\frac{f(\mathbf{y}+t\mathbf{d})-f(\mathbf{y})}{t}.$

For the class of locally Lipschitz functions,
$f^{\circ}(\mathbf{x};\mathbf{d})$ always exists for any
$\mathbf{x},\mathbf{d}\in\mathbb{R}^{n}$, and it is shown that
$f^{\circ}(\mathbf{x};\mathbf{d})=\max\\{\boldsymbol{\xi}^{T}\mathbf{d}\,\,:\,\,\boldsymbol{\xi}\in\partial
f(\mathbf{x})\\}$ [29].

The notion of _Pareto optimality_ is of crucial importance in multiobjective
optimization. A point $\mathbf{x}\in\mathbb{R}^{n}$ is called Pareto optimal
for problem (1) if there is not any $\mathbf{y}\in\mathbb{R}^{n}$ such that
$f_{i}(\mathbf{y})\leq f_{i}(\mathbf{x})$, for all $i=1,\ldots,p$, and
$f_{j}(\mathbf{y})<f_{j}(\mathbf{x})$, for some $j\in\\{1,\ldots,p\\}$ [2]. It
is recalled that for a point $\mathbf{x}\in\mathbb{R}^{n}$ to be a Pareto
optimal, it is necessary that [1]

$\mathbf{0}\in\texttt{conv}\left\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x})\right\\}.$

A point $\mathbf{x}$ which satisfies the above condition is called a _Clarke
substationary_ point.

In the current paper, the notion of weakly upper semismoothness is important
in designing an effective line search procedure.

######  Definition 2.1 ([30]).

The locally Lipschitz function $f:\mathbb{R}^{n}\to\mathbb{R}$ is called
weakly upper semismooth if for any $\mathbf{z},\mathbf{d}\in\mathbb{R}^{n}$
and sequences $\\{\boldsymbol{\xi}_{i}\\}_{i}\subset\mathbb{R}^{n}$ and
$\\{h_{i}\\}_{i}\subset\mathbb{R}_{+}$ satisfying $h_{i}\downarrow 0$ as
$i\to\infty$ and $\boldsymbol{\xi}_{i}\in\partial
f(\mathbf{z}+h_{i}\mathbf{d})$, one has

$\limsup_{i\to\infty}\boldsymbol{\xi}_{i}^{T}\mathbf{d}\geq\liminf_{i\to\infty}\frac{f(\mathbf{z}+h_{i}\mathbf{d})-f(\mathbf{z})}{h_{i}}.$

The class of weakly upper semismooth functions encompasses a wide range of
functions, including convex, concave, and functions of max- and min-type [30,
31, 32]. In fact, semismoothness is a popular assumption which is widely used
in the context of general nonsmooth nonconvex optimization [6, 10, 33].

##  3 Finding a descent direction

In this section, we develop an iterative procedure to find an efficient
descent direction for problem (1) at a given point
$\mathbf{x}\in\mathbb{R}^{n}$. A direction $\mathbf{d}\in\mathbb{R}^{n}$ is
called a descent direction for problem (1) at point
$\mathbf{x}\in\mathbb{R}^{n}$ if there exists $z>0$ such that

$f_{i}(\mathbf{x}+t\mathbf{d})-f_{i}(\mathbf{x})<0,\quad\text{for
all}\,\,t\in(0,z),\,\,\text{and}\,\,i=1,\ldots,p.$ (2)

The following lemma, which is essentially well-known, presents a theoretical
approach to find a descent direction. For the sake of completeness, we provide
a proof. Before it, let us denote the normal cone [34] of the convex set
$C\subset\mathbb{R}^{n}$ at a point $\mathbf{x}\in C$ by
$\mathcal{N}(C,\mathbf{x})$.

######  Lemma 3.1.

Assume $\boldsymbol{\xi}^{*}\neq\mathbf{0}$ solves the following minimization
problem:

$\min\big{\\{}\lVert\boldsymbol{\xi}\rVert\,\,:\,\,\boldsymbol{\xi}\in\text{\rm{conv}}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x})\\}\big{\\}}.$ (3)

Then the normalized direction

${\bar{\mathbf{d}}}:=-\boldsymbol{\xi}^{*}/\lVert\boldsymbol{\xi}^{*}\rVert,$
(4)

is a descent direction for multiobjective problem (1) at point $\mathbf{x}$.

###### Proof.

Compactness of $\texttt{conv}\\{\cup_{i=1}^{p}\partial f_{i}(\mathbf{x})\\}$
yields the existence of $\boldsymbol{\xi}^{*}$. Since $\boldsymbol{\xi}^{*}$
solves problem (3), we deduce [35]

$-\boldsymbol{\xi}^{*}\in\mathcal{N}\left(\texttt{conv}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x})\\},\boldsymbol{\xi}^{*}\right),$

yielding

$\boldsymbol{\xi}^{T}\boldsymbol{\xi}^{*}\geq\lVert\boldsymbol{\xi}^{*}\rVert^{2},\quad\text{for
all}\,\,\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x})\\},$

which is equivalent to

$-\boldsymbol{\xi}^{T}\frac{\boldsymbol{\xi}^{*}}{\lVert\boldsymbol{\xi}^{*}\rVert}\leq-\lVert\boldsymbol{\xi}^{*}\rVert<0,\quad\text{for
all}\,\,\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x})\\}.$

Consequently, $\boldsymbol{\xi}^{T}\mathbf{d}<0$ for all
$\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x})\\}$, which means $f^{\circ}_{i}(\mathbf{x};\mathbf{d})<0$,
for $i=1,\ldots,p$, and the assertion follows immediately from the definition
of $f^{\circ}(\mathbf{x};\mathbf{d})$. ∎

For some $i\in\\{1,\ldots,p\\}$, assume $\mathbf{x}\in\mathbb{R}^{n}$ is a
continuously differentiable point for the objective $f_{i}$, which is close to
a nonsmooth region. In this case, the subdifferential set $\partial
f_{i}(\mathbf{x})$ coincides with the singleton $\\{\nabla
f_{i}(\mathbf{x})\\}$. In fact, $\partial f_{i}(\mathbf{x})$ does not contain
any useful information of the nearby nonsmooth region, and consequently the
search direction ${\bar{\mathbf{d}}}$ given by (4) may not be an appropriate
descent direction. More precisely, when $\mathbf{x}$ is near a point where
$f_{i}$ is not differentiable, this search direction fulfills decrease
condition (2) only for very small values of $t$. In contrast,
$\partial_{\varepsilon}f_{i}(\mathbf{x})$ can capture some local information
of the nearby nonsmooth region, and one can stabilize our choice of search
direction by replacing $\partial f_{i}(\mathbf{x})$ by
$\partial_{\varepsilon}f_{i}(\mathbf{x})$. In this respect, the following
minimization problem is considered:

$\min\big{\\{}\lVert\boldsymbol{\xi}\rVert\,\,:\,\,\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial_{\varepsilon}f_{i}(\mathbf{x})\\}\big{\\}}.$
(5)

One major difficulty with problem (5) is that, for each objective $f_{i}$, we
need to know the entire subdifferential
$\partial_{\varepsilon}f_{i}(\mathbf{x})$, which is not an easy task in many
practical situations. In this regard, we develop an iterative algorithm to
efficiently approximate $\partial_{\varepsilon}f_{i}(\mathbf{x})$, for
$i=1,\ldots,p$.

For given $\varepsilon>0,m\in\mathbb{N}$ and $\mathbf{x}\in\mathbb{R}^{n}$,
let

$\mathbf{\mathcal{S}}_{\varepsilon}(\mathbf{x}):=\left[\begin{array}[]{cccc}\boldsymbol{\xi}_{1,1}&\boldsymbol{\xi}_{1,2}&\ldots&\boldsymbol{\xi}_{1,p}\\\
\boldsymbol{\xi}_{2,1}&\boldsymbol{\xi}_{2,2}&\ldots&\boldsymbol{\xi}_{2,p}\\\
\vdots&\vdots&\ddots&\vdots\\\
\boldsymbol{\xi}_{m,1}&\boldsymbol{\xi}_{m,2}&\ldots&\boldsymbol{\xi}_{m,p}\end{array}\right],$

be a collection of subgradients computed throughout iterations $1$ to $m$
having the property that

$\boldsymbol{\xi}_{j,i}\in\partial_{\varepsilon}f_{i}(\mathbf{x}),\quad\text{for
all}\,\,i=1,\ldots,p,\,\,\text{and}\,\,j=1,\ldots,m.$

In fact, the $i$-th column of $\mathbf{\mathcal{S}}_{\varepsilon}(\mathbf{x})$
stores a number of elements of $\partial_{\varepsilon}f_{i}(\mathbf{x})$. We
denote the $i$-th column of $\mathbf{\mathcal{S}}_{\varepsilon}(\mathbf{x})$
by $\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})$. Then, for each
$i\in\\{1,\ldots,p\\}$,

$\texttt{conv}\\{\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})\\}\subset\partial_{\varepsilon}f_{i}(\mathbf{x}),$

is an inner approximation of $\partial_{\varepsilon}f_{i}(\mathbf{x})$, and
consequently a practical variant of problem (5) can be given by

$\min\big{\\{}\lVert\boldsymbol{\xi}\rVert\,\,:\,\,\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})\\}\big{\\}}.$
(6)

Suppose $\boldsymbol{\xi}^{*}\neq\mathbf{0}$ solves problem (6), and define
the normalized direction ${\mathbf{d}}\in\mathbb{R}^{n}$ by

${\mathbf{d}}:=-\boldsymbol{\xi}^{*}/\lVert\boldsymbol{\xi}^{*}\rVert.$

If
$\texttt{conv}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})\\}$
is an adequate approximation of
$\texttt{conv}\\{\cup_{i=1}^{p}\partial_{\varepsilon}f_{i}(\mathbf{x})\\}$,
then one can find the step size $t>0$ satisfying

$f_{i}(\mathbf{x}+t{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta
t\lVert\boldsymbol{\xi}^{*}\rVert,\quad\text{for
all}\,\,i=1,\ldots,p,\quad\text{and}\quad t\geq\bar{t},$ (7)

in which $\beta\in(0,1)$ is the Armijo parameter which guarantees a sufficient
decrease in objectives $f_{i}$, and the lower bound
$\bar{t}\in(0,\varepsilon)$ excludes small values for the step size $t$. In
the case that condition (7) is met, we employ the search direction
$\mathbf{d}$ for taking a serious step, i.e., the current point $\mathbf{x}$
is updated by $\mathbf{x}^{+}:=\mathbf{x}+t\mathbf{d}.$ Otherwise, we infer
that our approximation of
$\texttt{conv}\\{\cup_{i=1}^{p}\partial_{\varepsilon}f_{i}(\mathbf{x})\\}$
must be improved. For this purpose, let

$\mathcal{I}:=\left\\{i\in\\{1,\ldots,p\\}\,\,:\,\,\nexists\,\,t\geq\bar{t}\,\,\,\,\text{\rm
s.t.}\,\,\,\,f_{i}(\mathbf{x}+t{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta
t\lVert\boldsymbol{\xi}^{*}\rVert\right\\}.$ (8)

Then, for each $i\in\mathcal{I}$, we need to improve our approximation of
$\partial_{\varepsilon}f_{i}(\mathbf{x})$. However, in practice, computing the
index set $\mathcal{I}$ is not possible because, for some $i$, checking the
condition

$\nexists\,\,t\geq\bar{t}\,\,\,\,\text{\rm
s.t.}\,\,\,\,f_{i}(\mathbf{x}+t{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta
t\lVert\boldsymbol{\xi}^{*}\rVert,$

is impractical. Therefore, in Algorithm 1, we propose a limited backtracking
line search to provide an estimation of the index set $\mathcal{I}$, namely
$\tilde{\mathcal{I}}$. Regarding this algorithm, some explanations are as
follows. The algorithm consists of two loops. The upper loop tests the
sufficient decrease condition (7) for a limited number of step sizes. If the
algorithm terminates within this loop, we conclude $\mathcal{I}=\emptyset$,
and hence $\tilde{\mathcal{I}}=\emptyset$ is returned. Moreover, one can
employ the obtained step size $t\geq\bar{t}>0$ to take a serious step.
Otherwise, the algorithm executes the lower loop, in which the step size
$\bar{t}>0$ is used to provide the nonempty estimation $\tilde{\mathcal{I}}$
of $\mathcal{I}$. Furthermore, the step size $t$ is set to zero as an
indication of taking a null step. Consequently, at termination of Algorithm 1,
one of the following cases occurs:

  * •

Case I. Algorithm 1 returns $\\{\tilde{\mathcal{I}},t\\}$ with
$\tilde{\mathcal{I}}=\emptyset$ and $t\geq\bar{t}>0$.

  * •

Case II. Algorithm 1 returns $\\{\tilde{\mathcal{I}},t\\}$ with
$\tilde{\mathcal{I}}\neq\emptyset$ and $t=0$.

Inputs: Objectives $\\{f_{i}\\}_{i=1}^{p}$, current point
$\mathbf{x}\in\mathbb{R}^{n}$, search direction
${\mathbf{d}}=-\boldsymbol{\xi}^{*}/\lVert\boldsymbol{\xi}^{*}\rVert$ with
$\boldsymbol{\xi}^{*}\neq\mathbf{0}$, Armijo coefficient $\beta\in(0,1)$, and
lower bound $\bar{t}\in(0,\varepsilon)$.

Parameters: Reduction factor $r\in(0,1)$, initial step size
$t_{0}\in(\bar{t},\infty)$, and number of backtracks $\tau\in\mathbb{N}$ with
$r^{\tau}t_{0}>\bar{t}>r^{\tau+1}t_{0}$.

Outputs: Step size $t\geq 0$ and index set
$\tilde{\mathcal{I}}\subset\\{1,\ldots,p\\}$.

Function: $\\{t,\tilde{\mathcal{I}}\\}$ = LBLS
($\\{f_{i}\\}_{i=1}^{p},\mathbf{x},\mathbf{d},\beta,\bar{t}$ )

1: Initialization: Set $\tilde{\mathcal{I}}:=\emptyset$ ;

2: for $t\in\\{t_{0},rt_{0},r^{2}t_{0},\ldots,r^{\tau}t_{0},\bar{t}\,\\}$ do

3: if $f_{i}(\mathbf{x}+t{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta
t\lVert\boldsymbol{\xi}^{*}\rVert$, for all $i=1,\ldots,p$, then

4: return {$t,\tilde{\mathcal{I}}$} and STOP ;

5: end if

6: end for

7: for $i\in\\{1,\ldots,p\\}$ do

8: if
$f_{i}(\mathbf{x}+\bar{t}{\mathbf{d}})-f_{i}(\mathbf{x})>-\beta\bar{t}\,\lVert\boldsymbol{\xi}^{*}\rVert$,
then

9: $\tilde{\mathcal{I}}:=\tilde{\mathcal{I}}\cup\\{i\\};$

10: end if

11: end for

12: Set $t:=0$ ;

13: return {$t,\tilde{\mathcal{I}}$} and STOP ;

End Function

Algorithm 1 Limited Backtracking Line Search (LBLS)

If Case I takes place, as already mentioned, we take a serious step by setting
$\mathbf{x}^{+}:=\mathbf{x}+t\mathbf{d}$. Thus, let us proceed with Case II.
In this case, we append the following row of subgradients

$[\boldsymbol{\xi}_{m+1,1},\boldsymbol{\xi}_{m+1,2},\ldots,\boldsymbol{\xi}_{m+1,p}],$

to the current collection $\mathbf{\mathcal{S}}_{\varepsilon}(\mathbf{x})$
such that

$\displaystyle\begin{split}&(\text{i})\,\,\boldsymbol{\xi}_{m+1,i}:=\boldsymbol{\xi}_{m,i},\,\,\text{for
all}\,\,i\notin\tilde{\mathcal{I}}.\\\
&(\text{ii})\,\,\boldsymbol{\xi}_{m+1,i}\in\partial_{\varepsilon}f_{i}(\mathbf{x})\,\,\text{and}\,\,\boldsymbol{\xi}_{m+1,i}\notin\texttt{conv}\\{\mathcal{S}_{\varepsilon,i}(\mathbf{x})\\},\,\,\text{for
all}\,\,i\in\tilde{\mathcal{I}}.\end{split}$ (9)

The first part of condition (9) states that for the objectives $f_{i}$ with
$i\notin\tilde{\mathcal{I}}$, we do not compute a new element of
$\partial_{\varepsilon}f_{i}(\mathbf{x})$. By the second part of this
condition, for each $i\in\tilde{\mathcal{I}}$, we compute a new element of
$\partial_{\varepsilon}f_{i}(\mathbf{x})$ which significantly improves our
approximation of $\partial_{\varepsilon}f_{i}(\mathbf{x})$. In other words, if
we update $\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})$ by

$\mathbf{\mathcal{S}}_{\varepsilon,i}^{+}(\mathbf{x}):=\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})\cup\\{\boldsymbol{\xi}_{m+1,i}\\},$

then part (ii) implies
$\texttt{conv}\\{\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})\\}\subsetneq\texttt{conv}\\{\mathbf{\mathcal{S}}_{\varepsilon,i}^{+}(\mathbf{x})\\}$.
Concerning part (ii) of condition (9), in the following lemma, we present a
sufficient condition for the new subgradient
$\boldsymbol{\xi}_{m+1,i}\in\partial_{\varepsilon}f_{i}(\mathbf{x})$ to
satisfy

$\boldsymbol{\xi}_{m+1,i}\notin\texttt{conv}\\{\mathcal{S}_{\varepsilon,i}(\mathbf{x})\\}.$

######  Lemma 3.2.

Assume $\bar{i}\in\tilde{\mathcal{I}}$ and
$\boldsymbol{\xi}^{*}\neq\mathbf{0}$ solves problem (6). For the direction
${\mathbf{d}}=-\boldsymbol{\xi}^{*}/\lVert\boldsymbol{\xi}^{*}\rVert$ and
scalar $c\in(0,1)$, suppose
$\boldsymbol{\xi}_{m+1,\bar{i}}\in\partial_{\varepsilon}f_{\bar{i}}(\mathbf{x})$
satisfies

$\boldsymbol{\xi}_{m+1,\bar{i}}^{T}{\mathbf{d}}\geq-c\lVert\boldsymbol{\xi}^{*}\rVert.$
(10)

Then
$\boldsymbol{\xi}_{m+1,\bar{i}}\notin\text{\rm{conv}}\\{\mathbf{\mathcal{S}}_{\varepsilon,\bar{i}}(\mathbf{x})\\}$.

###### Proof.

Since $\boldsymbol{\xi}^{*}$ is the optimal solution of problem (6), one can
conclude [35]

$-\boldsymbol{\xi}^{*}\in\mathcal{N}\left(\texttt{conv}\\{\cup_{i=1}^{p}\mathcal{S}_{\varepsilon,i}(\mathbf{x})\\},\boldsymbol{\xi}^{*}\right),$

which implies

$\boldsymbol{\xi}^{T}\boldsymbol{\xi}^{*}\geq\lVert\boldsymbol{\xi}^{*}\rVert^{2},\quad\text{for
all}\,\,\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\mathcal{S}_{\varepsilon,i}(\mathbf{x})\\}.$
(11)

Thus, if
$\boldsymbol{\xi}_{m+1,\bar{i}}\in\partial_{\varepsilon}f_{\bar{i}}(\mathbf{x})$
satisfies
$\boldsymbol{\xi}_{m+1,\bar{i}}^{T}{\mathbf{d}}\geq-c\lVert\boldsymbol{\xi}^{*}\rVert$,
inequality (11) yields
$\boldsymbol{\xi}_{m+1,\bar{i}}\notin\texttt{conv}\\{\cup_{i=1}^{p}\mathcal{S}_{\varepsilon,i}(\mathbf{x})\\}$,
and consequently

$\boldsymbol{\xi}_{m+1,\bar{i}}\notin\text{\rm{conv}}\\{\mathbf{\mathcal{S}}_{\varepsilon,\bar{i}}(\mathbf{x})\\}.$

∎

In Algorithm 2, we develop a new variant of Mifflin’s line search [6] to find
a new subgradient for objective $f_{i},\,i\in\tilde{\mathcal{I}}$, satisfying
condition (10). This algorithm consists of two conditional blocks. The first
one modifies the interval within which we look for an effective subgradient,
and the second one checks condition (10). It is also noted that, in this
algorithm, $t_{s}$ varies within the interval $(0,\varepsilon)$ and
$\boldsymbol{\xi}_{s}\in\partial f_{i}(\mathbf{x}+t_{s}\mathbf{d})$, for every
$s\geq 0$. Therefore, the output of this algorithm is an element of
$\partial_{\varepsilon}f_{i}(\mathbf{x})$.

Inputs: Objective $f_{i}$ with $i\in\tilde{\mathcal{I}}$, current point
$\mathbf{x}\in\mathbb{R}^{n}$, search direction
${\mathbf{d}}=-\boldsymbol{\xi}^{*}/\lVert\boldsymbol{\xi}^{*}\rVert$ with
$\boldsymbol{\xi}^{*}\neq\mathbf{0}$, and radius $\varepsilon\in(0,1)$.

Parameters: Armijo coefficient $\beta\in(0,1)$ and lower bound
$\bar{t}\in(0,\varepsilon)$ as set in Algorithm 1, scale parameter $c\in(0,1)$
with $\beta<c$, and reduction factor $\eta\in(0,0.5)$.

Output: A subgradient
$\boldsymbol{\xi}\in\partial_{\varepsilon}f_{i}(\mathbf{x})$.

Function: $\boldsymbol{\xi}$ = FES ($f_{i}$, $\mathbf{x}$, ${\mathbf{d}}$,
$\varepsilon$)

1: Initialization: Set $t_{0}:=\bar{t}$, compute
$\boldsymbol{\xi}_{0}\in\partial f_{i}(\mathbf{x}+t_{0}{\mathbf{d}})$, and set
$t^{l}_{0}:=0,t^{u}_{0}:=\varepsilon$, $s:=0$ ;

2: while true do

3: if
$f_{i}(\mathbf{x}+t_{s}{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta\,t_{s}\,\lVert\boldsymbol{\xi}^{*}\rVert,$
then

4: Set $t^{l}_{s+1}:=t_{s},\,\,\,\,t^{u}_{s+1}:=t^{u}_{s}$ ;

5: else

6: Set $t^{l}_{s+1}:=t^{l}_{s},\,\,\,\,t^{u}_{s+1}:=t_{s}$ ;

7: end if

8: if
$\boldsymbol{\xi}_{s}^{T}{\mathbf{d}}\geq-c\lVert\boldsymbol{\xi}^{*}\rVert,$
then

9: return $\boldsymbol{\xi}_{s}$ and STOP ;

10: end if

11: Choose
$t_{s+1}\in\left[t^{l}_{s+1}+\eta(t^{u}_{s+1}-t^{l}_{s+1}),\,t^{u}_{s+1}-\eta(t^{u}_{s+1}-t^{l}_{s+1})\right]$
;

12: Compute $\boldsymbol{\xi}_{s+1}\in\partial
f_{i}(\mathbf{x}+t_{s+1}{\mathbf{d}})$ ;

13: Set $s:=s+1$ ;

14: end while

End Function

Algorithm 2 Finding an Effective Subgradient (FES)

In the rest of this section, it is proved that Algorithm 2 terminates after
finitely many iterations. To this end, we provide some auxiliary results in
the next lemma.

######  Lemma 3.3.

If Algorithm 2 does not terminate ($s\to\infty$), then the followings hold:

  * (i)

We have $t_{s}\in\\{t_{s+1}^{l},t_{s+1}^{u}\\}$, for all $s\geq 0$. In
addition, for all $s\geq 1$,

$\displaystyle 0<t^{u}_{s+1}-t^{l}_{s+1}\leq(1-\eta)(t^{u}_{s}-t^{l}_{s}),$
(12) $\displaystyle 0\leq t^{l}_{s}\leq t^{l}_{s+1}<t^{u}_{s+1}\leq
t^{u}_{s}\leq\varepsilon.$ (13)

  * (ii)

There exists $t^{*}\in[0,\varepsilon]$ such that $t^{u}_{s}\downarrow t^{*}$,
$t^{l}_{s}\uparrow t^{*}$, and $t_{s}\to t^{*}$ as $s\to\infty$. Moreover

$t^{*}\in\mathbf{T}:=\\{t\,\,:\,\,f_{i}(\mathbf{x}+t{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta\,t\,\lVert\boldsymbol{\xi}^{*}\rVert\\}.$

  * (iii)

Suppose $\mathbf{S}:=\\{s\,\,:\,\,t^{u}_{s+1}=t_{s}\\}$. Then $\mathbf{S}$ is
an infinite set.

###### Proof.

(i) This part follows immediately from the construction of the algorithm.

(ii) It follows from (12) that the sequence $\\{t^{u}_{s}-t^{l}_{s}\\}_{s}$ is
decreasing and bounded by zero. Therefore, there exists $a\in\mathbb{R}$ such
that $(t^{u}_{s}-t^{l}_{s})\to a$ as $s\to\infty$. Letting $s$ approach
infinity in inequality (12), one can conclude $0\leq a\leq(1-\eta)a$. Since
$\eta\in(0,0.5)$, we deduce $a=0$. As $(t^{u}_{s}-t^{l}_{s})\to 0$ when
$s\to\infty$, inequalities (13) ensure the existence of
$t^{*}\in[0,\varepsilon]$ such that $t^{l}_{s}\uparrow
t^{*},t^{u}_{s}\downarrow t^{*}$ as $s\to\infty$. Moreover, since
$t_{s}\in\\{t^{l}_{s+1},t^{u}_{s+1}\\}$ for all $s\geq 0$, we have $t_{s}\to
t^{*}$ as $s\to\infty$. It remains to show that $t^{*}\in\mathbf{T}$. We have
$t^{l}_{s}\in\mathbf{T}$ for all $s\geq 0$, i.e.,

$f_{i}(\mathbf{x}+t^{l}_{s}{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta\,t^{l}_{s}\,\lVert\boldsymbol{\xi}^{*}\rVert,\qquad\text{for
all}\,\,s\geq 0.$

Now, continuity of $f_{i}$ along with the fact that $t^{l}_{s}\uparrow t^{*}$
as $s\to\infty$ yields

$f_{i}(\mathbf{x}+t^{*}{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta\,t^{*}\,\lVert\boldsymbol{\xi}^{*}\rVert,$

and therefore $t^{*}\in\mathbf{T}.$

(iii) First, we show $\mathbf{S}\neq\emptyset$. Assume by contradiction that
$\mathbf{S}=\emptyset$, which means

$f_{i}(\mathbf{x}+t_{s}{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta\,t_{s}\,\lVert\boldsymbol{\xi}^{*}\rVert,\qquad\text{for
all}\,\,s\geq 0.$ (14)

In particular, for $s=0$, we have $t_{0}=\bar{t}$, and thus

$f_{i}(\mathbf{x}+\bar{t}{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta\,\bar{t}\,\lVert\boldsymbol{\xi}^{*}\rVert,$

yielding $i\notin\tilde{\mathcal{I}}$ (see Algorithm 1), which is a
contradiction with the fact that $i\in\tilde{\mathcal{I}}$ in Algorithm 2.
Thus, $\mathbf{S}\neq\emptyset$. Next, we prove that $\mathbf{S}$ is indeed an
infinite set. Suppose by contradiction that $\mathbf{S}$ is finite. Then,
since $\mathbf{S}\neq\emptyset$ and $t^{u}_{s}\downarrow t^{*}$ as
$s\to\infty$, one can find $\bar{s}\in\mathbb{N}$ such that

$t^{u}_{s}=t^{*},\qquad\text{for all}\,\,s\geq\bar{s}\quad\text{and}\quad
t^{u}_{s}>t^{*},\qquad\text{for all}\,\,s<\bar{s}.$

Consequently, $t^{*}=t^{u}_{\bar{s}}=t_{\bar{s}-1}$, and therefore

$f_{i}(\mathbf{x}+t_{\bar{s}-1}{\mathbf{d}})-f_{i}(\mathbf{x})>-\beta\,t_{\bar{s}-1}\,\lVert\boldsymbol{\xi}^{*}\rVert,$

which gives $t^{*}\notin\mathbf{T}$, contradicting the fact
$t^{*}\in\mathbf{T}$. ∎

Now, we are in a position to state the principal convergence result for
Algorithm 2.

######  Theorem 3.4.

In Algorithm 2 suppose that objective function
$f_{i}:\mathbb{R}^{n}\to\mathbb{R}$ is weakly upper semismooth. Then this
algorithm terminates after a finite number of iterations.

###### Proof.

By indirect proof, assume that Algorithm 2 does not terminate. Let
$\mathbf{S}$ be as defined in part (iii) of Lemma 3.3. Then $\mathbf{S}$ is
infinite and

$f_{i}(\mathbf{x}+t_{s}{\mathbf{d}})-f_{i}(\mathbf{x})>-\beta
t_{s}\lVert\boldsymbol{\xi}^{*}\rVert,\quad\text{for all}\,\,s\in\mathbf{S}.$
(15)

By part (ii) of Lemma 3.3, we know that $t_{s}\to t^{*}$ as $s\to\infty$ with
$t^{*}\in\mathbf{T}$, i.e.,

$f_{i}(\mathbf{x}+t^{*}{\mathbf{d}})-f_{i}(\mathbf{x})\leq-\beta
t^{*}\lVert\boldsymbol{\xi}^{*}\rVert.$ (16)

Plugging the latter inequality into (15), we obtain

$f_{i}(\mathbf{x}+t_{s}{\mathbf{d}})-f_{i}(\mathbf{x}+t^{*}{\mathbf{d}})>-\beta\lVert\boldsymbol{\xi}^{*}\rVert(t_{s}-t^{*}),\quad\text{for
all}\,\,s\in\mathbf{S}.$ (17)

Define $\mathbf{z}:=\mathbf{x}+t^{*}{\mathbf{d}}$ and $h_{s}:=t_{s}-t^{*}>0$,
for all $s\in\mathbf{S}$. Then (17) can be rewritten as

$-\beta\lVert\boldsymbol{\xi}^{*}\rVert<\frac{f_{i}(\mathbf{z}+h_{s}{\mathbf{d}})-f_{i}(\mathbf{z})}{h_{s}},\quad\text{for
all}\,\,s\in\mathbf{S}.$ (18)

Since the objective $f_{i}$ is weakly upper semismooth, inequality (18)
implies

$-\beta\lVert\boldsymbol{\xi}^{*}\rVert\leq\liminf_{s\xrightarrow{\mathbf{S}}\infty}\frac{f_{i}(\mathbf{z}+h_{s}{\mathbf{d}})-f_{i}(\mathbf{z})}{h_{s}}\leq\limsup_{s\xrightarrow{\mathbf{S}}\infty}\boldsymbol{\xi}_{s}^{T}{\mathbf{d}}.$
(19)

On the other hand, since $s\to\infty$ in Algorithm 2, the stopping criterion
of this algorithm is never met, and hence

$\boldsymbol{\xi}_{s}^{T}{\mathbf{d}}<-c\lVert\boldsymbol{\xi}^{*}\rVert,\quad\text{for
all}\,\,s\in\mathbf{S},$

and consequently

$\limsup_{s\xrightarrow{\mathbf{S}}\infty}\boldsymbol{\xi}_{s}^{T}{\mathbf{d}}\leq-c\lVert\boldsymbol{\xi}^{*}\rVert<-\beta\lVert\boldsymbol{\xi}^{*}\rVert,$

which is a contradiction with (19). ∎

##  4 Finding a Clarke substationary point

The main goal of this section is to employ the provided tools in Section 3 in
order to locate a Clarke substationary point for multiobjective minimization
problem (1). For this purpose, the following concept of substationarity plays
a crucial role.

######  Definition 4.1 ([4, 5, 21]).

Suppose $\varepsilon>0$ and $\delta>0$ are given. A point
$\mathbf{x}\in\mathbb{R}^{n}$ is called a
$\big{(}\delta,\mathcal{S}_{\varepsilon}(\mathbf{x})\big{)}$-substationary
point for problem (1) if

$\min\big{\\{}\lVert\boldsymbol{\xi}\rVert\,\,:\,\,\boldsymbol{\xi}\in\text{\rm{conv}}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}_{\varepsilon,i}(\mathbf{x})\\}\big{\\}}\leq\delta.$

###  4.1 Computing a
$\big{(}\delta,\mathcal{S}_{\varepsilon}(\mathbf{x})\big{)}$-substationary
point

Based on the discussions provided in the previous section, a procedure for
finding a
$\big{(}\delta,\mathcal{S}_{\varepsilon}(\mathbf{x})\big{)}$-substationary
point is presented in Algorithm 3. Next, we aim to show that Algorithm 3
terminates after finite number of iterations. To this end, set

$\mathcal{K}:=\\{k\in\mathbb{N}_{0}\,\,:\,\,{t}_{k}>0\,\,\text{in Algorithm
\ref{Alg3}}\\}.$ (20)

Inputs: Objective functions $\\{f_{i}\\}_{i=1}^{p}$, starting point
$\mathbf{x}_{0}\in\mathbb{R}^{n}$, radius $\varepsilon\in(0,1)$,
substationarity tolerance $\delta>0$, and Armijo coefficient $\beta\in(0,1)$.

Parameter: Lower bound $\bar{t}\in(0,\varepsilon).$

Output: A $(\delta,\mathcal{S}_{\varepsilon}(\mathbf{x}))$-substationary point
$\mathbf{x}\in\mathbb{R}^{n}$.

Function: $\mathbf{x}$ = DS-SP ($\\{f_{i}\\}_{i=1}^{p}$,
$\mathbf{x}_{0},\varepsilon,\delta$, $\beta$)

1: Initialization: Compute $\boldsymbol{\xi}_{0,i}\in\partial
f_{i}(\mathbf{x}_{0})$, and set
$\mathcal{S}^{0}_{\varepsilon,i}(\mathbf{x}_{0}):=\\{\boldsymbol{\xi}_{0,i}\\}$,
for all $i\in\\{1,\ldots,p\\}$. Set $k:=0$ ;

2: while true do

3: Set
$\boldsymbol{\xi}^{*}_{k}:=\arg\min\big{\\{}\lVert\boldsymbol{\xi}\lVert\,\,:\,\,\boldsymbol{\xi}\in\text{\rm{conv}}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}^{k}_{\varepsilon,i}(\mathbf{x}_{k})\\}\big{\\}}$
;

4: if $\lVert\boldsymbol{\xi}^{*}_{k}\rVert\leq\delta,$ then

5: return $\mathbf{x}_{k}$ as a
$(\delta,\mathcal{S}_{\varepsilon}(\mathbf{x}_{k}))$-substationary point and
STOP ;

6: end if

7: Set
$\mathbf{d}_{k}:=-\boldsymbol{\xi}^{*}_{k}/\lVert\boldsymbol{\xi}^{*}_{k}\rVert$
;

8: Set $\\{t_{k},\tilde{\mathcal{I}_{k}}\\}:=$ LBLS
($\\{f_{i}\\}_{i=1}^{p},\mathbf{x}_{k},\mathbf{d}_{k},\beta,\bar{t}$ ) ;

9: if $t_{k}>0,$ then

10: Set $\mathbf{x}_{k+1}:=\mathbf{x}_{k}+t_{k}\mathbf{d}_{k}$ ;

11: Compute $\boldsymbol{\xi}_{k+1,i}\in\partial f_{i}(\mathbf{x}_{k+1})$, for
all $i\in\\{1,\ldots,p\\}$ ;

12: Set
$\mathcal{S}^{k+1}_{\varepsilon,i}(\mathbf{x}_{k+1}):=\\{\boldsymbol{\xi}_{k+1,i}\\}$,
for all $i\in\\{1,\ldots,p\\}$ ;

13: end if

14: if $t_{k}=0,$ then

15: Set
$\boldsymbol{\xi}_{k+1,i}:=\texttt{FES}(f_{i},\mathbf{x}_{k},\mathbf{d}_{k},\varepsilon)$,
for all $i\in\tilde{\mathcal{I}}_{k}$ ;

16: Set $\boldsymbol{\xi}_{k+1,i}:=\boldsymbol{\xi}_{k,i}$, for all
$i\notin\tilde{\mathcal{I}}_{k}$ ;

17: Set $\mathbf{x}_{k+1}:=\mathbf{x}_{k}$ ;

18: Set
$\mathcal{S}^{k+1}_{\varepsilon,i}(\mathbf{x}_{k+1}):=\mathcal{S}^{k}_{\varepsilon,i}(\mathbf{x}_{k})\cup\\{\boldsymbol{\xi}_{k+1,i}\\}$,
for all $i\in\\{1,\ldots,p\\}$ ;

19: end if

20: Set $k:=k+1$ ;

21: end while

End Function

Algorithm 3 Finding a
$\big{(}\delta,\mathcal{S}_{\varepsilon}(\mathbf{x})\big{)}$-Substationary
Point

We also make the following assumption about the multiobjective problem (1).

######  Assumption 4.1.

Objective functions $\\{f_{i}\\}_{i=1}^{p}$ are weakly upper semismooth.
Moreover, for some $j\in\\{1,\ldots,p\\}$, the sublevel set

$lev_{f_{j}(\mathbf{x}_{0})}(f_{j}):=\\{\mathbf{x}\in\mathbb{R}^{n}\,\,:\,\,f_{j}(\mathbf{x})\leq
f_{j}({\mathbf{x}_{0}})\\},$

is bounded.

We start with the following lemma.

######  Lemma 4.1.

Suppose that Assumption 4.1 holds. If Algorithm 3 does not terminate, then the
index set $\mathcal{K}$ is finite.

###### Proof.

Since $f_{j}$ is continuous, $lev_{f_{j}(\mathbf{x}_{0})}(f_{j})$ is closed,
and hence, by assumption, $lev_{f_{j}(\mathbf{x}_{0})}(f_{j})$ is compact.
Therefore

$f^{*}_{j}:=\min\\{f_{j}(\mathbf{x})\,\,:\,\,\mathbf{x}\in\mathbb{R}^{n}\\}>-\infty.$
(21)

By indirect proof, suppose that $\mathcal{K}$ is infinite. As Algorithm 3 does
not terminate, we conclude

$\lVert\boldsymbol{\xi}^{*}_{k}\rVert>\delta,\quad\text{for all}\,\,k.$ (22)

For each $k\in\mathcal{K}$, we have ${t}_{k}>0$ and thus

$f_{i}(\mathbf{x}_{k+1})-f_{i}(\mathbf{x}_{k})\leq-\beta\,t_{k}\,\lVert\boldsymbol{\xi}^{*}_{k}\rVert,\quad\text{for
all}\,\,i=1,\ldots,p\,\,\,\text{and}\,\,\,k\in\mathcal{K}.$

By construction of Algorithm 1, we have $t_{k}\geq\bar{t}>0$, for all
$k\in\mathcal{K}$. Thus, in virtue of (22), the latter inequality implies

$f_{i}(\mathbf{x}_{k+1})-f_{i}(\mathbf{x}_{k})\leq-\beta\,\bar{t}\,\delta,\quad\text{for
all}\,\,i=1,\ldots,p\,\,\,\text{and}\,\,\,k\in\mathcal{K}.$

Furthermore, for any $k\in\mathbb{N}_{0}\setminus\mathcal{K}$, we have
${t_{k}}=0$ and hence

$f_{i}(\mathbf{x}_{k+1})=f_{i}(\mathbf{x}_{k}),\quad\text{for
all}\,\,i=1,\ldots,p\,\,\,\text{and}\,\,\,k\in\mathbb{N}_{0}\setminus\mathcal{K}.$

In particular, for $i=j$, one can write

$\displaystyle
f_{j}(\mathbf{x}_{k+1})-f_{j}(\mathbf{x}_{k})\leq-\beta\,\bar{t}\,\delta,\quad\text{for
all}\,\,k\in\mathcal{K},$ (23) $\displaystyle
f_{j}(\mathbf{x}_{k+1})=f_{j}(\mathbf{x}_{k}),\quad\text{for
all}\,\,k\in\mathbb{N}_{0}\setminus\mathcal{K}.$ (24)

Using (23) and (24) inductively, we arrive at

$f_{j}(\mathbf{x}_{k+1})\leq
f_{j}(\mathbf{x}_{0})-\sum_{\begin{subarray}{c}q\in\mathcal{K}\\\ q\leq
k+1\end{subarray}}\beta\,\bar{t}\,\delta.$ (25)

As $\mathcal{K}$ is infinite, $\sum_{\begin{subarray}{c}q\in\mathcal{K}\\\
q\leq k+1\end{subarray}}\beta\,\bar{t}\,\delta\to\infty$ as $k\to\infty$.
Thus, (25) gives $f_{j}(\mathbf{x}_{k})\to-\infty$ as $k\to\infty$, which
contradicts (21). ∎

Our main result about Algorithm 3 is given in the next theorem.

######  Theorem 4.2.

Suppose that Assumption 4.1 holds. Then Algorithm 3 terminates in finite
number of iterations.

###### Proof.

By indirect proof, suppose that Algorithm 3 does not terminate, i.e.,
$k\to\infty$. Therefore

$\lVert\boldsymbol{\xi}_{k}^{*}\lVert>\delta,\quad\text{for all}\,\,k.$ (26)

Let $\mathcal{K}$ be as defined in (20). Then, in view of Lemma 4.1,
$\mathcal{K}$ is a finite set. In the case that $\mathcal{K}\neq\emptyset$,
set

$\bar{k}:=\max\\{k\,\,:\,\,k\in\mathcal{K}\\}.$

Otherwise, we set $\bar{k}:=0$. Define
$\bar{\mathbf{x}}:=\mathbf{x}_{\bar{k}+1}$. Then, for any $k>\bar{k}$, we have
${t}_{k}=0$, and consequently

$\bar{\mathbf{x}}=\mathbf{x}_{k},\quad\text{for all}\,\,k>\bar{k}.$

Furthermore, for any $k>\bar{k}$, the following collection of new subgradients
is computed:

$\boldsymbol{\chi}_{k}:=\\{\boldsymbol{\xi}_{k+1,i}\,\,:\,\,\boldsymbol{\xi}_{k+1,i}=\texttt{FES}(f_{i},\bar{\mathbf{x}},\mathbf{d}_{k},\varepsilon)\in\partial_{\varepsilon}f_{i}(\bar{\mathbf{x}}),\,\,i\in\tilde{\mathcal{I}}_{k}\\}.$

By construction of Algorithm 2, each element of $\boldsymbol{\chi}_{k}$
satisfies

$\boldsymbol{\xi}_{k+1,i}^{T}{\mathbf{d}}_{k}\geq-c\lVert\boldsymbol{\xi}^{*}_{k}\lVert.$

Since
${\mathbf{d}}_{k}=-\boldsymbol{\xi}^{*}_{k}/\lVert\boldsymbol{\xi}^{*}_{k}\lVert$,
the above inequality is represented as

$\boldsymbol{\xi}_{k+1,i}^{T}\boldsymbol{\xi}^{*}_{k}\leq
c\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}.$ (27)

For all $k>\bar{k}$, we have $\tilde{\mathcal{I}}_{k}\neq\emptyset$. Let
$i_{k}\in\tilde{\mathcal{I}}_{k}$ be arbitrary. Then,
$\boldsymbol{\xi}_{k+1,i_{k}}\in\boldsymbol{\chi}_{k}$, and hence

$\boldsymbol{\xi}_{k+1,i_{k}}\in\mathcal{S}^{k+1}_{\varepsilon,i_{k}}(\mathbf{x}_{k+1})=\mathcal{S}^{k+1}_{\varepsilon,i_{k}}(\bar{\mathbf{x}})\subset\partial_{\varepsilon}f_{i_{k}}(\bar{\mathbf{x}})\subset\texttt{conv}\\{\cup_{i=1}^{p}\partial_{\varepsilon}f_{i}(\bar{\mathbf{x}})\\}.$
(28)

It is also noted that, for all $k>\bar{k}$,

$\boldsymbol{\xi}^{*}_{k}\in\texttt{conv}\\{\cup_{i=1}^{p}\mathcal{S}^{k}_{\varepsilon,i}({\mathbf{x}_{k}})\\}=\texttt{conv}\\{\cup_{i=1}^{p}\mathcal{S}^{k}_{\varepsilon,i}(\bar{\mathbf{x}})\\}\subset\texttt{conv}\\{\cup_{i=1}^{p}\partial_{\varepsilon}f_{i}(\bar{\mathbf{x}})\\}.$
(29)

Boundedness of
$\texttt{conv}\\{\cup_{i=1}^{p}\partial_{\varepsilon}f_{i}(\bar{\mathbf{x}})\\}$
implies

$C:=\sup\\{\lVert\boldsymbol{\xi}\rVert\,\,:\,\,\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial_{\varepsilon}f_{i}(\bar{\mathbf{x}})\\}\\}<\infty.$

Set $\bar{C}:=\max\\{C,\delta\\}$. Therefore, it follows from (28) and (29)
that

$\lVert\boldsymbol{\xi}_{k+1,i_{k}}-\boldsymbol{\xi}^{*}_{k}\rVert\leq
2\bar{C},\quad\text{for all}\,\,k>\bar{k}.$ (30)

Now, for any $t\in(0,1)$ and $k>\bar{k}$, one can write

$\displaystyle\lVert\boldsymbol{\xi}^{*}_{k+1}\rVert^{2}$
$\displaystyle\leq\lVert
t\boldsymbol{\xi}_{k+1,i_{k}}+(1-t)\boldsymbol{\xi}^{*}_{k}\lVert^{2}$
$\displaystyle=t^{2}\lVert\boldsymbol{\xi}_{k+1,i_{k}}-\boldsymbol{\xi}^{*}_{k}\lVert^{2}+2t(\boldsymbol{\xi}^{*}_{k})^{T}(\boldsymbol{\xi}_{k+1,i_{k}}-\boldsymbol{\xi}^{*}_{k})+\lVert\boldsymbol{\xi}^{*}_{k}\rVert^{2}.$
(31)

In view of (27) and (30), one can continue (4.1) as

$\displaystyle\lVert\boldsymbol{\xi}^{*}_{k+1}\rVert^{2}$ $\displaystyle\leq
4t^{2}\bar{C}^{2}+2tc\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}-2t\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}+\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}$
$\displaystyle=4t^{2}\bar{C}^{2}+\left(1-2t(1-c)\right)\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}$
$\displaystyle=:\varphi(t),$ (32)

for all $t\in(0,1)$. It is easy to see that

$t^{*}:=2(1-c)\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}/8\bar{C}^{2}\in(0,1),$

is the minimum point of $\varphi(t)$ and

$\varphi(t^{*})=\left(1-\frac{(1-c)^{2}\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}}{4\bar{C}^{2}}\right)\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}.$

In virtue of (26), the above equality gives

$\varphi(t^{*})\leq\left(1-\frac{(1-c)^{2}\delta^{2}}{4\bar{C}^{2}}\right)\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2}.$
(33)

Since $\delta\leq\bar{C}$ and $c\in(0,1)$, we get

$\alpha:=1-\frac{(1-c)^{2}\delta^{2}}{4\bar{C}^{2}}\in(0,1).$

Next, it follows from (32) and (33) that

$0\leq\lVert\boldsymbol{\xi}^{*}_{k+1}\rVert^{2}\leq\varphi(t^{*})\leq\alpha\lVert\boldsymbol{\xi}^{*}_{k}\lVert^{2},\quad\text{for
all}\,\,k>\bar{k},$ (34)

yielding the sequence
$\\{\lVert\boldsymbol{\xi}^{*}_{k}\rVert^{2}\\}_{k>\bar{k}}$ is bounded from
below and decreasing. Thus, it converges. Let
$\lVert\boldsymbol{\xi}^{*}_{k}\rVert^{2}\to A$ as $k\to\infty$. Letting
$k\to\infty$ in inequality (32), we obtain $0\leq A\leq\alpha A$. Since
$\alpha\in(0,1)$, we deduce $A=0$. Consequently,
$\lVert\boldsymbol{\xi}^{*}_{k}\rVert^{2}\to 0$ as $k\to\infty$, which
contradicts (26).

∎

###  4.2 Computing a Clarke substationary point

In this section, for the given decreasing sequences $\delta_{\nu}\downarrow 0$
and $\varepsilon_{\nu}\downarrow 0$, we seek a Clarke substationary point of
problem (1) by generating a sequence of
$\big{(}\delta_{\nu},\mathcal{S}_{\varepsilon_{\nu}}(\cdot)\big{)}$-substationary
points. Such a sequence is simply generated by Algorithm 4.

Inputs: Objective functions $\\{f_{i}\\}_{i=1}^{p}$, starting point
$\mathbf{x}_{0}\in\mathbb{R}^{n}$, positive sequences $\delta_{\nu}\downarrow
0$ and $\varepsilon_{\nu}\downarrow 0$, and optimality tolerance $\rho>0$.

Parameter: Armijo coefficient $\beta>0$.

Output: A point $\mathbf{x}\in\mathbb{R}^{n}$ as an approximation of a Clarke
substationary point.

1: Initialization: Set $\nu:=0$ ;

2: while true do

3: Set $\mathbf{x}_{\nu+1}:=\,\texttt{DS-
SP}\,(\\{f_{i}\\}_{i=1}^{p},\mathbf{x}_{\nu},\varepsilon_{\nu},\delta_{\nu},\beta)$
;

4: if $\delta_{\nu}<\rho,$ and $\varepsilon_{\nu}<\rho,$ then

5: return $\mathbf{x}_{\nu}$ as an approximation of a Clarke substationary
point and Stop ;

6: end if

7: Set $\nu:=\nu+1$ ;

8: end while

Algorithm 4 Computation of a Clarke stationary point

To study the asymptotic behavior of Algorithm 4, one may assume $\rho=0$. In
this case, the algorithm generates the infinite sequence
$\\{\mathbf{x}_{\nu}\\}_{\nu}$. In the next theorem, it is proved that any
accumulation point of the sequence $\\{\mathbf{x}_{\nu}\\}_{\nu}$ is Clarke
substationary for problem (1).

######  Theorem 4.3.

Suppose that Assumption 4.1 holds. If $\rho=0$ in Algorithm 4, then any
accumulation point of the sequence $\\{\mathbf{x}_{\nu}\\}_{\nu}$ generated by
this algorithm is Clarke substationary for problem (1).

###### Proof.

During each iteration $\nu\geq 0$ of Algorithm 4, the DS-SP algorithm is run
and terminated after finite number of iterations. For any $\nu\geq 0$, suppose
that this algorithm is terminated at iteration $k_{\nu}\in\mathbb{N}_{0}$.
Thus

$\mathbf{x}_{\nu+1}=\mathbf{x}_{k_{\nu}}=\,\texttt{DS-
SP}\,(\\{f_{i}\\}_{i=1}^{p},\mathbf{x}_{\nu},\varepsilon_{\nu},\delta_{\nu},\beta),$

and therefore

$\min\big{\\{}\lVert\boldsymbol{\xi}\lVert\,\,:\,\,\boldsymbol{\xi}\in\text{\rm{conv}}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}^{k_{\nu}}_{\varepsilon_{\nu},i}(\mathbf{x}_{\nu+1})\\}\big{\\}}\leq\delta_{\nu},\quad\text{for
all}\,\,\nu\geq 0.$ (35)

We note that $\mathbf{x}_{\nu}\in lev_{f_{j}(\mathbf{x}_{0})}(f_{j})$, for all
$\nu\geq 0$, and by Assumption 4.1, the sub-level set
$lev_{f_{j}(\mathbf{x}_{0})}(f_{j})$ is bounded. Consequently, the sequence
$\\{\mathbf{x}_{\nu}\\}_{\nu}$ has at least one accumulation point, say
$\mathbf{x}^{*}$, i.e., there exists $\mathcal{V}\subset\mathbb{N}_{0}$ such
that $\mathbf{x}_{\nu}\xrightarrow{\mathcal{V}}\mathbf{x}^{*}$. This fact
along with (35) implies

$\min\big{\\{}\lVert\boldsymbol{\xi}\lVert\,\,:\,\,\boldsymbol{\xi}\in\text{\rm{conv}}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}^{k_{\nu}}_{\varepsilon_{\nu},i}(\mathbf{x}_{\nu+1})\\}\big{\\}}\leq\delta_{\nu},\quad\text{for
all}\,\,\nu\in\mathcal{V}.$ (36)

Let $\kappa>0$ be arbitrary. As $\delta_{\nu}\downarrow 0$ as $\nu\to\infty$,
one can find $\bar{\nu}\in\mathcal{V}$ sufficiently large such that
$\delta_{\nu}<\kappa$, for all $\nu\geq\bar{\nu}$. Let
$\mathcal{V}^{\prime}:=\\{\nu\in\mathcal{V}\,:\,\nu\geq\bar{\nu}\\}\subset\mathcal{V}$.
Then, it immediately follows from (36) that, for all
$\nu\in\mathcal{V}^{\prime}$,

$\lVert\boldsymbol{\xi}^{*}_{\nu}\rVert:=\min\big{\\{}\lVert\boldsymbol{\xi}\lVert\,\,:\,\,\boldsymbol{\xi}\in\text{\rm{conv}}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}^{k_{\nu}}_{\varepsilon_{\nu},i}(\mathbf{x}_{\nu+1})\\}\big{\\}}\leq\kappa.$
(37)

In addition, since

$\boldsymbol{\xi}^{*}_{\nu}\in\text{\rm{conv}}\\{\cup_{i=1}^{p}\mathbf{\mathcal{S}}^{k_{\nu}}_{\varepsilon_{\nu},i}(\mathbf{x}_{\nu+1})\\}\subset\text{\rm{conv}}\\{\cup_{i=1}^{p}\partial_{\varepsilon_{\nu}}f_{i}(\mathbf{x}_{\nu+1})\\},\quad\text{for
all}\,\,\nu\in\mathcal{V}^{\prime},$

there exist

$\boldsymbol{\xi}^{*}_{\nu,i}\in\partial_{\varepsilon_{\nu}}f_{i}(\mathbf{x}_{\nu+1}),\,\,\lambda_{\nu,i}\in\mathbb{R},\quad
i=1,\ldots,p,$

satisfying

$\boldsymbol{\xi}^{*}_{\nu}=\sum_{i=1}^{p}\lambda_{\nu,i}\boldsymbol{\xi}^{*}_{\nu,i},\quad\sum_{i=1}^{p}\lambda_{\nu,i}=1,\quad\lambda_{\nu,i}\geq
0,\,\,\,i=1,\ldots,p,$ (38)

for all $\nu\in\mathcal{V}^{\prime}$. The fact that
$\mathbf{x}_{\nu}\xrightarrow{\mathcal{V}^{\prime}}\mathbf{x}^{*}$ along with
uniformly compactness of the map $\partial_{\cdot}f_{i}(\cdot)$ implies that
the sequences
$\\{\boldsymbol{\xi}^{*}_{\nu,i}\\}_{\nu\in\mathcal{V}^{\prime}}$, for
$i=1,\ldots,p$, are bounded. Now, boundedness of the sequences
$\\{\lambda_{\nu,i}\\}_{\nu\in\mathcal{V}^{\prime}}$, for $i=1,\ldots,p$,
implies the existence of the infinite set
$\mathcal{V}^{\prime\prime}\subset\mathcal{V}^{\prime}$ such that

$\boldsymbol{\xi}^{*}_{\nu,i}\xrightarrow{\mathcal{V}^{\prime\prime}}\boldsymbol{\xi}^{*}_{i},\quad\lambda_{\nu,i}\xrightarrow{\mathcal{V}^{\prime\prime}}\lambda_{i},\quad\lambda_{i}\geq
0,\quad\sum_{i=1}^{p}\lambda_{i}=1,\quad\text{for all}\,\,i=1,\ldots,p.$ (39)

Furthermore, as the map $\partial_{\cdot}f_{i}(\cdot)$ is upper
semicontinuous, one can deduce

$\boldsymbol{\xi}^{*}_{i}\in\partial f_{i}(\mathbf{x}^{*}),\quad\text{for
all}\,\,i=1,\ldots,p.$ (40)

Eventually, in view of (38), (39) and (40), we conclude

$\boldsymbol{\xi}^{*}_{\nu}=\sum_{i=1}^{p}\lambda_{\nu,i}\boldsymbol{\xi}^{*}_{\nu,i}\xrightarrow{\mathcal{V}^{\prime\prime}}\sum_{i=1}^{p}\lambda_{i}\boldsymbol{\xi}^{*}_{i}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x}^{*})\\}.$ (41)

Consequently, (37) and (41) yield

$\min\big{\\{}\lVert\boldsymbol{\xi}\rVert\,\,:\,\,\boldsymbol{\xi}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x}^{*})\\}\big{\\}}\leq\kappa.$

Since $\kappa>0$ was arbitrary, we conclude
$\mathbf{0}\in\texttt{conv}\\{\cup_{i=1}^{p}\partial
f_{i}(\mathbf{x}^{*})\\}$. ∎

######  Remark 4.1.

Concerning Remark 4.1 in [21], our proof of Theorem 4.3 shows the convergence
of Algorithm 4 of Gebken and Peitz [21] to a Clarke substationary point when
$\delta_{\nu}\downarrow 0,\varepsilon_{\nu}\downarrow 0$ as $\nu\to\infty$.

##  5 Numerical experiments

In this section, we assess the practical performance of the proposed multi-
objective optimization method. To this end, first we analyze the typical
behavior of the proposed method. Next, we apply our method to approximate the
Pareto fronts of five multiobjective test problems. Then, we seek a sparse
solution to an underdetermined linear system by considering a nonsmooth bi-
objective minimization problem. Eventually, we apply the proposed method to a
collection of nonsmooth convex and nonconvex test problems, and some
comparative results with a state-of-the-art algorithm are provided. Throughout
this section, the proposed method is called SUMOPT (SUbgradient method for
Multiobjective OPTimization). We have implemented the method in Matlab
(R2017b) on a machine with CPU 3.5 GHz and 8 GB of RAM.

Regarding the parameters for SUMOPT, our choices are as follows.

  * •

Algorithm 1. In the LBLS algorithm, the reduction factor $r\in(0,1)$ and the
initial step size $t_{0}\in(\bar{t},\infty)$ are set to be $0.5$ and 2.0,
respectively. In addition, for the number of backtrack steps
$\tau\in\mathbb{N}$, we choose
$\tau:=\lceil{\frac{\ln{\bar{t}}-\ln{t_{0}}}{\ln{r}}-1}\rceil$.

  * •

Algorithm 2. In the FES algorithm, we set $c:=0.01$. Moreover, since
$\eta\in(0,0.5)$, in Line 11 of this algorithm, one can update $t_{s}$ by

$t_{s+1}:=0.5(t^{u}+t^{l}).$

  * •

Algorithm 3. In this algorithm, the lower bound $\bar{t}\in(0,\varepsilon)$ is
set to $\varepsilon/10$. Also, the quadratic subproblem in this algorithm is
solved using the _quadprog_ solver in Matlab.

  * •

Algorithm 4. We work with the Armijo coefficient $\beta:=10^{-6}$, and the
positive sequences $\delta_{\nu+1}:=\gamma\delta_{\nu}$,
$\varepsilon_{\nu+1}:=\gamma\varepsilon_{\nu}$ with
$\gamma=\delta_{0}=\varepsilon_{0}:=0.1$. Moreover, the starting point
$\mathbf{x}_{0}\in\mathbb{R}^{n}$ and optimality tolerance $\rho>0$ will be
specified in each experiment.

The parameters mentioned earlier were utilized in all the subsequent
experiments, except for subsection 5.1. In this particular case, we ran the
SUMOPT method using $\gamma=0.5$, $\delta_{0}=0.3,\bar{t}:={\varepsilon}/2$,
and $t_{0}=0.25$ to increase the number of iterations and obtain a more
comprehensive understanding of the method’s performance.

In the following experiments, we consider a collection of convex and nonconvex
multiobjective test problems with $n=2$ and $p\in\\{2,3,4,5\\}$, which are
presented in Table 1. In this table, “P” stands for the problem number, and
the symbols “+” and “-” denote convex and nonconvex objectives, respectively.
For more details on the considered objective functions, see [26, 36]. The
remaining test problems are provided and described in the corresponding
subsections.

Table 1: List of test problems

P |  | $f_{i}$ |  | Convexity |
---|---|---|---|---|---
1 |  | [Crescent, LQ] |  | [-, +] |
2 |  | [Mifflin2, Crescent] |  | [-, -] |
3 |  | [Crescent, QL] |  | [-, +] |
4 |  | [CB3, LQ] |  | [+, +] |
5 |  | [CB3, Mifflin1] |  | [+, +] |
6 |  | [Mifflin2, Mifflin1] |  | [-, +] |
7 |  | [CB3, QL] |  | [+, +] |
8 |  | [Mifflin2, DEM] |  | [-, +] |
9 |  | [Mifflin2, LQ] |  | [-, +] |
10 |  | [CB3, DEM] |  | [+, +] |
11 |  | [DEM, QL, Mifflin1] |  | [+, +, +] |
12 |  | [Mifflin2, Crescent, Mifflin1] |  | [-, -, +] |
13 |  | [DEM, QL, Mifflin1, CB3] |  | [+, +, +, +] |
14 |  | [Mifflin2, Crescent, DEM, Mifflin1] |  | [-, -, +, +] |
15 |  | [Mifflin2, Crescent, DEM, Mifflin1, QL] |  | [-, -, +, +, +] |

Throughout this section, we also consider a state-of-the-art algorithm, namely
the Multiobjective Proximal Bundle method (MPB) [22, 37] to provide some
comparative results. For this purpose, we used the Fortran code of the MPB
method from [37] using the default parameters.

###  5.1 Tracking the algorithm

In this experiment, we analyze the typical behavior of the proposed SUMOPT
algorithm using the test problem P1. Indeed, we consider the following bi-
objective test problem:

$\min\,\,\left(f_{1}(\mathbf{x}),f_{2}(\mathbf{x})\right)\quad\text{s.t.}\quad\mathbf{x}\in\mathbb{R}^{2},$
(42)

in which

$\displaystyle
f_{1}(\mathbf{x}):=\max\\{x_{1}^{2}+(x_{2}-1)^{2}+x_{2}-1,-x_{1}^{2}-(x_{2}-1)^{2}+x_{2}+1\\}\quad\text{(Crescent),}$
$\displaystyle
f_{2}(\mathbf{x}):=\max\\{-x_{1}-x_{2},-x_{1}-x_{2}+x_{1}^{2}+x_{2}^{2}-1\\}\quad\text{(LQ).}$

It is evident that $f_{1}(\mathbf{x})$ is a nonconvex function, while
$f_{2}(\mathbf{x})$ is convex, and the nonsmooth regions of these functions
are given by

$\mathcal{M}_{1}:=\\{\mathbf{x}\in\mathbb{R}^{2}\,:\,x_{1}^{2}+(x_{2}-1)^{2}=1\\}\quad\text{and}\quad\mathcal{M}_{2}:=\\{\mathbf{x}\in\mathbb{R}^{2}\,:\,x_{1}^{2}+x_{2}^{2}=1\\},$

respectively. Table 2 reports the observed results using the starting point
$\mathbf{x}_{0}^{T}=[-0.6,0.2]$ and the optimality tolerance $\rho=5\times
10^{-3}$. Each block of the table displays the results of Algorithm 3, which
are generated by Algorithm 4 at iteration $\nu\geq 0$.

Table 2: Numerical results of SUMOPT applied to problem (42).

$\nu=0\qquad\delta_{0}=0.3\qquad\varepsilon_{0}=0.1$
---
$k$ |  | $\lVert\boldsymbol{\xi}^{*}_{k}\rVert$ |  | $\mathbf{d}_{k}^{T}$ |  | $\tilde{I}_{k}$ |  | $\mathbf{x}_{k+1}^{T}$ |  | $f_{1}(\mathbf{x}_{k+1})$ |  | $f_{2}(\mathbf{x}_{k+1})$
0 |  | 1.3416 |  | (0.8944, 0.4472) |  | $\\{\texttt{1}\\}$ |  | (-0.6000, 0.2000) |  | 0.2000 |  | 0.4000
1 |  | 0.3494 |  | (0.8598, -0.5104) |  | $\varnothing$ |  | (-0.3850, 0.0723) |  | 0.0811 |  | 0.3126
2 |  | 1.1508 |  | (0.6691, 0.7431) |  | $\\{\texttt{1}\\}$ |  | (-0.3850, 0.0723) |  | 0.0811 |  | 0.3126
3 |  | 0.3925 |  | (0.9268, -0.3755) |  | $\varnothing$ |  | (-0.1533, -0.0214) |  | 0.0454 |  | 0.1748
4 |  | 1.0871 |  | (0.2820, 0.9593) |  | $\varnothing$ |  | (-0.1445, 0.0084) |  | 0.0124 |  | 0.1360
5 |  | 1.0246 |  | (0.2820, 0.9593) |  | $\\{\texttt{1}\\}$ |  | (-0.1445, 0.0084) |  | 0.0124 |  | 0.1360
6 |  | 0.1478 |  | – |  | – |  | – |  | – |  | –
$\nu=1\qquad\delta_{1}=0.15\qquad\varepsilon_{1}=0.05$
0 |  | 1.0246 |  | (0.2820, 0.9593) |  | $\\{\texttt{1}\\}$ |  | (-0.1445, 0.0084) |  | 0.0124 |  | 0.01360
1 |  | 0.1460 |  | – |  | – |  | – |  | – |  | –
$\nu=2\qquad\delta_{2}=0.075\qquad\varepsilon_{2}=0.025$
0 |  | 1.0246 |  | (0.2820, 0.9593) |  | $\varnothing$ |  | (-0.1439, 0.0103) |  | 0.0104 |  | 0.1335
1 |  | 1.0207 |  | (0.2820, 0.9593) |  | $\\{\texttt{1}\\}$ |  | (-0.1439, 0.0103) |  | 0.0104 |  | 0.1335
2 |  | 0.1448 |  | (0.9897, -0.1430) |  | $\varnothing$ |  | (-0.0224, -0.0075) |  | 0.0079 |  | 0.0277
3 |  | 1.0158 |  | (0.0398, 0.9992) |  | $\varnothing$ |  | (-0.0199, 0.0002) |  | 0.0004 |  | 0.0196
4 |  | 0.7161 |  | (0.9678, -0.2516) |  | $\\{\texttt{1}\\}$ |  | (-0.0199, 0.0002) |  | 0.0004 |  | 0.0196
5 |  | 0.0017 |  | – |  | – |  | – |  | – |  | –
$\nu=3\qquad\delta_{3}=0.0375\qquad\varepsilon_{3}=0.0125$
0 |  | 0.7161 |  | (0.9678, -0.2516) |  | $\varnothing$ |  | (-0.0189, 0.0000) |  | 0.0003 |  | 0.0189
1 |  | 1.0006 |  | (0.0379, 0.9992) |  | $\\{\texttt{1}\\}$ |  | (-0.0189, 0.0000) |  | 0.0003 |  | 0.0189
2 |  | 0.0190 |  | – |  | – |  | – |  | – |  | –
$\nu=4\qquad\delta_{4}=0.01875\qquad\varepsilon_{4}=0.00625$
0 |  | 1.0006 |  | (0.0379, 0.9992) |  | $\\{\texttt{1}\\}$ |  | (-0.0189, -0.0000) |  | 0.0003 |  | 0.0189
1 |  | 0.0190 |  | (0.9998, -0.0189) |  | $\varnothing$ |  | (-0.0033, -0.0002) |  | 0.0002 |  | 0.0036
2 |  | 1.0005 |  | (0.0067, 0.9999) |  | $\\{\texttt{1}\\}$ |  | (-0.0033, -0.0002) |  | 0.0002 |  | 0.0036
3 |  | 0.0033 |  | – |  | – |  | – |  | – |  | –
$\nu=5\qquad\delta_{5}=0.009375\qquad\varepsilon_{5}=0.003125$
0 |  | 1.0005 |  | (0.0067, 0.9999) |  | $\varnothing$ |  | (-0.0033, -0.0000) |  | 0.0000 |  | 0.0033
1 |  | 1.0000 |  | (0.0067, 0.9999) |  | $\\{\texttt{1}\\}$ |  | (-0.0033, -0.0000) |  | 0.0000 |  | 0.0033
2 |  | 0.0033 |  | – |  | – |  | – |  | – |  | –
$\nu=6\qquad\delta_{6}=0.0046875\qquad\varepsilon_{6}=0.0015625$
Termination

For $\nu=0$, since $\mathbf{x}_{0}\in\mathcal{M}_{1}$, we observe
$\tilde{I}_{0}=\\{1\\}$, which means the singleton set
$\mathcal{S}_{0.1,1}^{0}$ is not an adequate approximation of
$\partial_{0.1}f_{1}(\mathbf{x}_{0})$ to make a reduction in objective
function $f_{1}$. Consequently, the method admits a null step to enrich
$\mathcal{S}_{0.1,1}^{0}$ by using the FES algorithm. Note that, as
$2\notin\tilde{I}_{0}$, the LBLS algorithm has already obtained a descent step
for $f_{2}$ at $\mathbf{x}_{0}$, and hence we do not improve our approximation
of $\partial_{0.1}f_{2}(\mathbf{x}_{0})$ at iteration $k=0$. Once a new
effective element of $\partial_{0.1}f_{1}(\mathbf{x}_{0})$ is appended to
$\mathcal{S}_{0.1,1}^{0}$ to form $\mathcal{S}_{0.1,1}^{1}$, one can see
$\tilde{I}_{1}$ becomes empty, which means a serious step occurs at $k=1$.
Moreover, since $\mathcal{S}_{0.1,1}^{0}\subsetneq\mathcal{S}_{0.1,1}^{1}$, we
have
$\lVert\boldsymbol{\xi}^{*}_{1}\rVert<\lVert\boldsymbol{\xi}^{*}_{0}\rVert$.
In general, after each serious step, $\lVert\boldsymbol{\xi}^{*}_{k}\rVert$
becomes smaller. However, it should be noted that the sequence
$\\{\lVert\boldsymbol{\xi}^{*}_{k}\rVert\\}_{k}$ is not monotonically
decreasing, yet, as shown in the proof of Theorem 4.3, it has a decreasing
subsequence. Eventually, at iteration $k=6$,
$\lVert\boldsymbol{\xi}^{*}_{6}\rVert$ drops below $\delta_{0}=0.3$, and thus
Algorithm 3 terminates.

Regarding $\nu=1$, we do not see any serious step. In fact, after a null step,
the method recognizes that the current point $\mathbf{x}_{1}$ is indeed a
$\big{(}0.15,\mathcal{S}_{0.05}(\mathbf{x}_{1})\big{)}$-substationary point
(see Definition 4.1).

For $\nu=2$, Algorithm 3 is executed with starting point
$\mathbf{x}_{0}^{T}=(-0.1445,0.0084)$, at which the search direction
$\mathbf{d}_{0}^{T}=(0.2820,0.9593)$ is obtained. Although no serious step is
taken along this direction at $\nu=1$, for $\nu=2$ we observe
$\tilde{I}_{0}=\varnothing$, indicating that a serious step is taken. This is
because as $\nu$ increases, the lower bound $\bar{t}=\varepsilon_{\nu}/{2}$ in
the LBLS algorithm decreases, resulting in more backtrack steps. Consequently,
a serious step occurs. After this serious step, the current point
$\mathbf{x}_{1}$ approaches the nonsmooth region $\mathcal{M}_{1}$. Hence,
taking a null step allows the FES algorithm to capture a subgradient from the
opposite side of $\mathcal{M}_{1}$, which significantly modifies the search
direction $\mathbf{d}^{T}_{1}$. In fact, the modified search direction
$\mathbf{d}^{T}_{2}$ is almost perpendicular to $\mathbf{d}^{T}_{1}$,
resulting in an effective serious step. The method continues in this manner
until, at $\nu=6$, we observe $\delta_{6}$ and $\varepsilon_{6}$ become
smaller than the optimality tolerance $\rho=0.5\times 10^{-3}$, and the SUMOPT
algorithm terminates.

###  5.2 Approximating the Pareto front

In addition to convergence, the efficiency of a method in approximating the
Pareto front is an important property. Thus, the main goal of this experiment
is to employ the multi-start technique [38, 39] in order to generate an
approximation of the Pareto front. To this end, we consider problems P1-P5
given in Table 1. To generate different Pareto points using the SUMOPT and MPB
algorithms, we independently and uniformly generate $3\times 10^{2}$ random
starting points from the following two-dimensional box:

$\mathbf{B}:=\\{(x_{1},x_{2})\in\mathbb{R}^{2}\,\,:\,\,0\leq x_{1},x_{2}\leq
2\\}.$

Moreover, the optimality tolerance $\rho$ in the SUMOPT method is set to
$10^{-4}$. Figure 1 depicts the objective space of problems P1-P5 along with
the obtained approximations of the Pareto fronts by the SUMOPT (top plots) and
MPB (bottom plots) methods.

Figure 1: Objective space and obtained approximations of the Pareto fronts by
the SUMOPT (top plots) and MPB (bottom plots) methods, for problems P1-P5.

Once the Pareto set of a bi-objective optimization problem is connected, the
Hole Absolute Size (HAS) and Hole Relative Size (HRS) metrics [40] are
effective tools for assessing the uniformity of the points along the
approximate Pareto front. Assume that $\mathcal{Y}$ is a Pareto set
approximation whose elements have been sorted in ascending order based on the
first objective function. For $j=1,\ldots,\lvert\mathcal{Y}\rvert-1$, take
$d_{j}:=\lVert\mathbf{y}_{j}-\mathbf{y}_{j+1}\rVert$, which is the Euclidean
distance between two consecutive solutions
$\mathbf{y}_{j},\mathbf{y}_{j+1}\in\mathcal{Y}$. Then, if $\bar{D}$ is the
mean of all $d_{j}$ for $j=1,\ldots,\lvert\mathcal{Y}\rvert-1$, the indicators
HAS and HRS are defined by

$\displaystyle\textbf{HAS}(\mathcal{Y}):=\max\\{d_{j}\,\,:\,\,j=1,\ldots,\lvert\mathcal{Y}\rvert-1\\},$
$\displaystyle\textbf{HRS}(\mathcal{Y}):=\frac{1}{\bar{D}}\max\\{d_{j}\,\,:\,\,j=1,\ldots,\lvert\mathcal{Y}\rvert-1\\}.$

The HAS measures the absolute size of the largest gap between consecutive
solutions, while the HRS provides a relative measure of this gap compared to
the average gap size. Obviously, for these indicators, a lower value is
desirable.

Table 3 presents the metrics values and the cardinality of the approximate
Pareto fronts, $\lvert\mathcal{Y}\rvert$, for the obtained approximations
addressed in Figure 1. By comparing the values of HAS and HRS between the two
methods, one can conclude that SUMOPT offers a more uniform approximation of
the Pareto fronts in the considered set of test problems. In addition, it can
be seen that for problem P2 and the MPB method, sixty-two of the starting
points did not result in Pareto points. The same happened for problem P3 and
the SUMOPT method for seven of the starting points.

Table 3: Metrics HAS and HRS along with $\lvert\mathcal{Y}\rvert$, for
problems P1-P5.

|  |  | SUMOPT |  |  | MPB |
---|---|---|---|---|---|---|---
P |  |  | HAS |  | HRS |  | $\lvert\mathcal{Y}\rvert$ |  |  | HAS |  | HRS |  | $\lvert\mathcal{Y}\rvert$ |
1 |  |  | 0.0952 |  | 13.6061 |  | 300 |  |  | 0.1403 |  | 23.3442 |  | 300 |
2 |  |  | 0.1379 |  | 29.0251 |  | 300 |  |  | 0.1846 |  | 31.7277 |  | 238 |
3 |  |  | 1.5664 |  | 9.2399 |  | 293 |  |  | 2.4199 |  | 14.7110 |  | 300 |
4 |  |  | 0.0544 |  | 11.6060 |  | 300 |  |  | 0.1334 |  | 34.2662 |  | 300 |
5 |  |  | 0.6107 |  | 8.7263 |  | 300 |  |  | 1.6905 |  | 23.5538 |  | 300 |

Our next example presents an interesting case. It involves a smooth bi-
objective problem where the objective space is not
$\mathbb{R}^{p}_{\geqq}$-convex [2], and the set of Pareto front is not
connected. In such a scenario, it is well-known that the Weighted Sum (WS)
method [2, 15], as a popular scalarization technique, falls short of
generating a suitable approximation of the Pareto front. Let us design the
following bi-objective test problem:

$\min\,\,\left(f_{1}(x),f_{2}(x)\right)\quad\text{s.t.}\quad x\in\mathbb{R},$
(43)

where $f_{1},f_{2}:\mathbb{R}\to\mathbb{R}$ are given by

$\displaystyle
f_{1}(x):=\big{(}1+0.1\sin(8x)\big{)}\big{(}\cos(x)\cos(0.6)-\sin(-0.6)\sin(x)\big{)},$
$\displaystyle
f_{2}(x):=\big{(}1+0.1\sin(8x)\big{)}\big{(}\cos(x)\sin(-0.6)+\sin(x)\cos(0.6)\big{)}.$

Since $f_{1}$ and $f_{2}$ are periodic functions with periodicity $T=2\pi$,
the objective space of problem (43) can be obtained, e.g., on the interval
$I=[0,2\pi]$, see Figure 2 (top-left). Let $G(x):=\nabla f_{1}(x)\nabla
f_{2}(x)$. Then, it is easy to verify that the set of substationary points of
problem (43) on the interval $I$ is given by ${S}:=\\{x\in I\,\,:\,\,G(x)\leq
0\\}$. The top-middle plot of Figure 2 depicts the graph of $G(x)$ together
with its nonpositive parts, $G^{-}(x)$, associated with the set $S$.
Furthermore, the set of substationary points along with the set of Pareto
points in the objective space have been shown in the right-top and left-bottom
plots of the figure, respectively. To achieve an approximation of the Pareto
front, we applied the SUMOPT method using a uniform grid of $10^{3}$ starting
points from the interval $[0,2\pi]$, and the optimality tolerance
$\rho=10^{-3}$. One can see from the bottom-middle plot of Figure 2 that
SUMOPT has successfully generated the entire Pareto elements.

Figure 2: Top: details of the considered test problem (43). Bottom: Pareto
points and obtained approximation of Pareto elements by SUMOPT and WS methods.

For problem (43), we consider the weighted sum problem, which is formulated as
follows:

$\min_{x\in\mathbb{R}}\lambda f_{1}(x)+(1-\lambda)f_{2}(x),$ (44)

where the weight $\lambda$ varies in the interval $[0,1]$. To solve this
problem for various weights, we selected $\lambda$ from a uniform grid of
$10^{3}$ points within the interval $[0,1]$. Next, for each $\lambda$, problem
(43) was solved using the _patternsearch_ solver in Matlab. The obtained
results are depicted in the bottom-right plot of Figure 2. It is noteworthy
that the WS technique fails to generate a significant portion of the Pareto
front.

###  5.3 Sparse solution to an underdetermined linear system

In this experiment, we consider the following linear system

$\mathbf{A}\mathbf{x}=\mathbf{b},$ (45)

in which, $\mathbf{A}\in\mathbb{R}^{m\times n}$ with $m<n$, and
$\mathbf{b}\in\mathbb{R}^{m}$. This system is referred to as underdetermined
and typically has an infinite number of solutions. Such systems appear in
various practical problems, including data fitting, machine learning,
denoising, and image processing. The main goal is to find an approximate
solution to linear system (45) that achieves the highest level of sparsity.
The classic approach to formulate this problem is as follows:

$\displaystyle\min\,\,\lVert\mathbf{A}\mathbf{x}-\mathbf{b}\rVert_{2}^{2}+\lambda\lVert\mathbf{x}\rVert_{1},\quad\text{s.t.}\quad\mathbf{x}\in\mathbb{R}^{n},$
(46)

where $\lVert\mathbf{x}\rVert_{1}$ denotes the $l_{1}$ norm of the vector
$\mathbf{x}\in\mathbb{R}^{n}$, and $\lambda>0$ is a positive parameter.
However, choosing an appropriate value for the parameter $\lambda$ is a
delicate and challenging issue in effectively solving this problem [41]. In
this regard, one alternative interesting approach is to formulate the problem
by the following bi-objective optimization problem:

$\displaystyle\min\,\,(\lVert\mathbf{x}\rVert_{1},\lVert\mathbf{A}\mathbf{x}-\mathbf{b}\rVert_{2}^{2})\quad\text{s.t.}\quad\mathbf{x}\in\mathbb{R}^{n}.$
(47)

This problem is convex, but nonsmooth due to the use of $l_{1}$ norm. To
create an instance of problem (47), we randomly generated the matrix
$\mathbf{A}\in\mathbb{R}^{50\times 100}$ and vector
$\mathbf{b}\in\mathbb{R}^{50}$. The SUMOPT was then applied to this problem
using five randomly generated starting points from a unit box, and the
optimality tolerance $\rho=10^{-2}$. For each run, the obtained results have
been reported in Table 4.

Table 4: Results of the SUMOPT method on problem (47), for five randomly
generated starting points.

Run number |  | $\lVert\mathbf{x}^{*}\rVert_{1}$ |  | $\lVert\mathbf{A}\mathbf{x}^{*}-\mathbf{b}\rVert_{2}^{2}$ |  | #Zero |  | Time (s) |  | #Fun |  | #Sub |
---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 |  | 1.7632 |  | 1.2932 |  | 76 |  | 4.7154 |  | 14788 |  | 748 |
2 |  | 0.4483 |  | 5.4864 |  | 91 |  | 14.6328 |  | 50002 |  | 2189 |
3 |  | 0.7555 |  | 2.8640 |  | 92 |  | 8.8045 |  | 30008 |  | 1367 |
4 |  | 1.1938 |  | 1.9667 |  | 82 |  | 3.9671 |  | 13542 |  | 680 |
5 |  | 2.1281 |  | 1.0854 |  | 85 |  | 4.7224 |  | 16182 |  | 797 |

In this table, $\mathbf{x}^{*}\in\mathbb{R}^{100}$ is the obtained approximate
Pareto point. Additionally, #Fun and #Sub denote the total number of function
and subgradient evaluations during the minimization process, respectively.
Furthermore, #Zero indicates the number of components in $\mathbf{x}^{*}$ with
absolute values less than $10^{-3}$. Indeed, we define a zero component as any
having an absolute value less than $10^{-3}$. For all runs, SUMOPT reached the
desired accuracy within less than fifteen seconds of CPU time. According to
the table, we can see that for the second and third runs, we obtained
solutions with fewer non-zero elements, but this comes at a cost of increased
residual $\lVert\mathbf{A}\mathbf{x}^{*}-\mathbf{b}\rVert_{2}^{2}$. Table 4
also highlights an interesting trade-off: the solution with the smallest
$l_{1}$ norm has the largest residual, and the solution with the smallest
residual has the largest $l_{1}$ norm. Regarding the evaluation of functions
and subgradients, we can observe that the number of function evaluations
exceeds the number of subgradient evaluations. This is because the SUMOPT
algorithm uses a limited backtracking line search to eliminate objectives that
do not require computation of a new subgradient after a null step.

###  5.4 A comparison with the MPB method

This experiment aims to compare the efficiency of the proposed SUMOPT
algorithm with the MPB method in terms of the number of function and
subgradient evaluations. To this end, we applied these methods to a collection
of fifteen nonsmooth convex and nonconvex multiobjective test problems
presented in Table 1. Regarding the MPB method, we used the Fortran code from
[37] with the default parameters. For the SUMOPT method, we worked with the
aforementioned parameters with the optimality tolerance $\rho=10^{-3}$.

Concerning the set of starting points, following [21], we considered a
$13\times 13$ uniform grid of the region
$[-3,3]\times[-3,3]\subset\mathbb{R}^{2}$ as follows:

$\mathbf{S}:=\\{(x_{1},x_{2})\in\mathbb{R}^{2}\,\,:\,\,x_{1}=-3+0.5i,\,\,x_{2}=-3+0.5j,\,\,i,j=0,1,\ldots,12\\}.$

Indeed, we ran each problem 169 times using the set of starting points
$\mathbf{S}$, and the total number of evaluations is reported in Table 5. It
should be noted that MPB and SUMOPT have their own different first-order
optimality condition. Yet, by this choice of parameters, both methods
successfully reached the desired accuracy in all runs, and the observed
results are distinguishable enough so that one can highlight the prominent
feature of each.

In Table 5, for SUMOPT, “Iter” represents the total number of iterations of
Algorithm 3 during the runs, while for MPB, it denotes the total number of
serious steps performed across the runs. Additionally, #Fun and #Sub indicate
the total number of function and subgradient evaluations throughout the runs.
For each objective, MPB evaluates the function value and one arbitrary
subgradient at each trial point. Consequently, the number of function and
subgradient evaluations is always the same in this method.

As observed from the table, for all test problems, MPB requires fewer function
evaluations than our proposed method. This can be attributed to the fact that
SUMOPT performs a limited backtracking line search to exclude objectives for
which new subgradient computations are not required after a null step.
However, SUMOPT is more efficient than MPB in terms of subgradient
evaluations, as MPB computes a new subgradient for all objectives after each
null step, whereas SUMOPT avoids unnecessary subgradient computations.
Regarding iterations, the same pattern as for the number of subgradient
evaluations occurs. In summary, SUMOPT achieves better subgradient evaluation
efficiency compared to MPB, though MPB requires fewer function evaluations
through its avoidance of backtracking steps.

Table 5: Results of the SUMOPT and MPB methods on the set of problems
considered in Table 1.

|  |  SUMOPT |  |  | MPB |
---|---|---|---|---|---|---
P |  | Iter |  | #Fun |  | #Sub |  |  | Iter |  | #Fun |  | #Sub |
1 |  | 1085 |  | 20560 |  | 2704 |  |  | 1723 |  | 4426 |  | 4426 |
2 |  | 878 |  | 13902 |  | 2367 |  |  | 1706 |  | 4528 |  | 4528 |
3 |  | 711 |  | 11030 |  | 1831 |  |  | 1894 |  | 4454 |  | 4454 |
4 |  | 660 |  | 7842 |  | 1973 |  |  | 1046 |  | 2634 |  | 2634 |
5 |  | 1724 |  | 30006 |  | 3761 |  |  | 2950 |  | 7332 |  | 7332 |
6 |  | 1086 |  | 18854 |  | 2655 |  |  | 2609 |  | 6842 |  | 6842 |
7 |  | 1205 |  | 20300 |  | 2777 |  |  | 1844 |  | 4068 |  | 4068 |
8 |  | 764 |  | 13146 |  | 2188 |  |  | 869 |  | 2118 |  | 2118 |
9 |  | 1247 |  | 24144 |  | 3050 |  |  | 1795 |  | 4352 |  | 4352 |
10 |  | 944 |  | 17314 |  | 2394 |  |  | 916 |  | 2278 |  | 2278 |
11 |  | 442 |  | 9050 |  | 2380 |  |  | 1078 |  | 3972 |  | 3972 |
12 |  | 853 |  | 19299 |  | 3516 |  |  | 2241 |  | 11733 |  | 11733 |
13 |  | 485 |  | 14004 |  | 3386 |  |  | 1006 |  | 4904 |  | 4904 |
14 |  | 654 |  | 18234 |  | 4289 |  |  | 1253 |  | 9088 |  | 9088 |
15 |  | 575 |  | 22721 |  | 4878 |  |  | 819 |  | 5070 |  | 5070 |

##  6 Conclusion

We have developed a descent subgradient method for finding a Clarke
substationary point of a general multiobjective optimization problem, whose
objectives are locally Lipschitz. As opposed to the multiobjective proximal
bundle method proposed by Mäkelä et al., we do not turn the main problem into
a single-objective one by using an improvement function. Instead, in order to
find an effective descent direction for the main problem, the behavior of each
individual objective is separately studied. This framework allows us to make
use of a limited backtracking line search, which helps us determine the
objectives for which improving the current approximation of the
$\varepsilon$-subdifferential set is unnecessary. This strategy reduced the
number of subgradient evaluations during the optimization process.

Regarding the enrichment process of the current collection of the subgradient
information, the proposed method shares a common feature with bundle-type
methods, in that the $\varepsilon$-subdifferential sets are sequentially
improved. In this regard, we presented a new variant of Mifflin’s line search
that provides a non-redundant subgradient, which significantly modifies the
current working set associated with each objective in the improvement process.
We proved that such a subgradient element can be found through a finite
procedure by establishing a relationship between the limited backtracking line
search and the proposed Mifflin’s line search. We studied the global
convergence of the method to a Clarke substationary point, which strengthened
the results derived in the method proposed by Gebken and Peitz.

In our numerical tests, we assessed the typical behavior of our algorithm,
SUMOPT, and observed the efficiency of the proposed limited backtracking line
search in excluding objective functions from the improvement process.
Meanwhile, the proposed Mifflin’s line search showed its ability in modifying
the inefficient search directions, after each null step. Next, we showed the
efficiency of our method in generating a uniform approximation of the Pareto
fronts. For finding a sparse solution of an underdetermined linear system, our
method revealed a trade-off between the sparsity of the solution and the value
of the remainder along the Pareto front of the problem. Eventually, our
comparative results based on fifteen nonsmooth test problems showed that
SUMOPT is a big rival to bundle-like methods, especially in case the
computation of subgradients is a time-consuming process.

## References

  * [1] K. Miettinen, _Nonlinear Multiobjective Optimization_. New York: Springer, 1999.
  * [2] M. Ehrgott, _Multicriteria optimization_. New York: Springer, 2005.
  * [3] Y. Sawaragi, H. Nakayama, and T. Tanino, _Theory of Multiobjective Optimization_. Waltham: Academic Press, 1985.
  * [4] A. M. Bagirov, L. Jin, N. Karmitsa, A. Al Nuaimat, and N. Sultanova, “A subgradient method for nonconvex nonsmooth optimization,” _J. Optim. Theory Appl._ , vol. 157, pp. 416–435, 2013.
  * [5] A. M. Bagirov and A. N. Ganjehlou, “A quasisecant method for minimizing nonsmooth functions,” _Optim. Methods Softw._ , vol. 25, no. (1), pp. 3–18, 2010.
  * [6] K. C. Kiwiel, _Methods of Descent for Nondifferentiable Optimization_. Berlin: Springer-Verlag, 1985.
  * [7] W. Hare, C. Sagastizábal, and M. Solodov, “A proximal bundle method for nonsmooth nonconvex functions with inexact information,” _Comput. Optim. Appl._ , vol. 63, no. (1), pp. 1–28, 2016.
  * [8] M. Maleknia and M. Shamsi, “A new method based on the proximal bundle idea and gradient sampling technique for minimizing nonsmooth convex functions,” _Comput. Optim. Appl._ , vol. 77, pp. 379–409, 2020.
  * [9] ——, “A quasi-Newton proximal bundle method using gradient sampling technique for minimizing nonsmooth convex functions,” _Optim. Methods Softw_ , vol. 37, no. 4, pp. 1415–1446, 2022.
  * [10] H. Schramm and J. Zowe, “A version of the bundle idea for minimizing a nonsmooth function: conceptual idea, convergence analysis, numerical results,” _SIAM J. Optim._ , vol. 2, no. (1), pp. 105–122, 1992.
  * [11] J. V. Burke, A. S. Lewis, and M. L. Overton, “A robust gradient sampling algorithm for nonsmooth, nonconvex optimization,” _SIAM J. Optim._ , vol. 15, no. (3), pp. 751–779, 2005.
  * [12] K. C. Kiwiel, “Convergence of the gradient sampling algorithm for nonsmooth nonconvex optimization,” _SIAM J. Optim._ , vol. 18, no. (2), pp. 379–388, 2007.
  * [13] F. E. Curtis and M. L. Overton, “A sequential quadratic programming algorithm for nonconvex, nonsmooth constrained optimization,” _SIAM J. Optim._ , vol. 22, no. (2), pp. 474–500, 2012.
  * [14] M. Maleknia and M. Shamsi, “A gradient sampling method based on ideal direction for solving nonsmooth optimization problems,” _J. Optim. Theory Appl._ , vol. 187, pp. 181–204, 2020.
  * [15] A. M. Geoffrion, “Proper efficiency and the theory of vector maximization,” _J. Optim. Theory Appl._ , vol. 22, pp. 618–630, 1968.
  * [16] M. Gravel, I. M. Martel, R. Madeau, W. Price, and R. Tremblay, “A multicriterion view of optimal ressource allocation in job-shop production,” _European J. Oper. Res._ , vol. 61, pp. 230–244, 1992.
  * [17] J. Fliege, L. M. G. Drummond, and B. Svaiter, “Newton’s method for multiobjective optimization,” _SIAM J. Optim_ , vol. 20, no. 2, pp. 602–626, 2009.
  * [18] J. Y. B. Cruz, “A subgradient method for vector optimization problems,” _SIAM J. Optim._ , vol. 23, no. 4, p. 2169–2182, 2013.
  * [19] J. X. D. C. Neto, G. J. P. D. Silva, O. P. Ferreira, and J. O. Lopes, “A subgradient method for multiobjective optimization,” _Comput Optim Appl_ , vol. 54, p. 461–472, 2013.
  * [20] O. Montonen, N. Karmitsa, and M. M. Mäkelä, “Multiple subgradient descent bundle method for convex nonsmooth multiobjective optimization,” _Comput Optim Appl_ , vol. 67, no. 1, pp. 139–158, 2018.
  * [21] B. Gebken and S. Peitz, “An efficient descent method for locally lipschitz multiobjective optimization problems,” _J. Optim. Theory Appl._ , vol. 188, no. 3, p. 696–723, 2021.
  * [22] M. M. Mäkelä, N. Karmitsa, and O. Wilppu, “Multiobjective proximal bundle method for nonsmooth optimization. TUCS technical report No 1120,” Turku Centre for Computer Science, Turku, Tech. Rep., 2014.
  * [23] S. Qu, C. Liu, M. Goh, Y. Li, and Y. Ji, “Nonsmooth multiobjective programming with quasi-Newton methods,” _European J. Oper. Res._ , vol. 235, no. 3, pp. 503–510, 2014.
  * [24] N. H. Monjezi and S. Nobakhtian, “An inexact multiple proximal bundle algorithm for nonsmooth nonconvex multiobjective optimization problems,” _Ann. Oper. Res._ , vol. 311, no. 2, pp. 1123–1154, 2022.
  * [25] ——, “A proximal bundle-based algorithm for nonsmooth constrained multiobjective optimization problems with inexact data,” _Numerical Algorithms_ , vol. 89, no. 2, pp. 637–674, 2022.
  * [26] A. M. Bagirov, N. Karmitsa, and M. M. Mäkelä, _Introduction to Nonsmooth Optimization_. Springer International Publishing, 2014.
  * [27] L. Lukšan and J. Vlček, “Test problems for nonsmooth unconstrained and linearly constrained optimization,” Institute of Computer Science, Academy of Sciences of the Czech Republic, Prague, Tech. Rep., 2000.
  * [28] L. C. Evans and R. F. Gariepy, _Measure Theory and Fine Properties of Functions, Revised Edition_. Boca Raton: CRC Press, 1992.
  * [29] F. H. Clarke, _Optimization and Nonsmooth Analysis_. Philadelphia: SIAM, 1990.
  * [30] R. Mifflin, “Semismooth and semiconvex functions in constrained optimization,” _SIAM J. Optim._ , vol. 15, no. 6, pp. 959–972, 1977.
  * [31] ——, “An algorithm for constrained optimization with semismooth functions,” _Mathematics of Operations Research_ , vol. 2, no. (2), pp. 191–207, 1977.
  * [32] A. M. Bagirov, N. Karmitsa, and S. Taheri, _Partitional Clustering via Nonsmooth Optimization_. Switzerland: Springer Cham, 2020.
  * [33] J. Outrata, M. Kocvara, and J. Zowe, _Nonsmooth Approach to Optimization Problems with Equilibrium Constraints_. Dordrecht: Springer, 1998.
  * [34] B. S. Mordukhovich, _Variational Analysis and Generalized Differentiation I_. Berlin: Springer, 2006.
  * [35] R. T. Rockafellar and R. J.-B. Wets, _Variational Analysis_. Berlin: Springer, 2004.
  * [36] M. M. Mäkelä and P. Neittaanmäki, _Nonsmooth Optimization: Analysis and Algorithms with Applications to Optimal Control_. Singapore: World Scientific Publishing Co., 1992.
  * [37] M. M. Mäkelä, “Multiobjective proximal bundle method for nonconvex nonsmooth optimization: Fortran Subroutine MPBNGC 2.0,” Depart. Math. Inf. Technol. Ser. B. Sci. Comput. B 13, 2003, Tech. Rep., 2003.
  * [38] J. Fliege, L. M. G. Drummond, and B. F. Svaiter, “Newton’s method for multiobjective optimization,” _SIAM J. Optim._ , vol. 20, no. 2, pp. 602–626, 2009.
  * [39] M. A. T. Ansary and G. Panda, “A globally convergent SQCQP method for multiobjective optimization problems,” _SIAM J. Optim._ , vol. 31, no. 1, pp. 91–113, 2021.
  * [40] C. Audet, J. Bigeon, D. Cartier, S. L. Digabel, and L. Salomon, “Performance indicators in multiobjective optimization,” _European J. Oper. Res._ , vol. 292, no. 2, pp. 397–422, 2021.
  * [41] A. Beck, _First-Order Methods in Optimization_. Philadelphia: SIAM, 2017.