# On the Treatment of Optimization Problems with L1 Penalty Terms via
Multiobjective Continuation

Katharina Bieker, Bennet Gebken and Sebastian Peitz K. Bieker and B. Gebken
are with the Department of Mathematics, Paderborn University, Germany.
E-mail<EMAIL_ADDRESS>
S. Peitz is with the Department of Computer Science, Paderborn University,
Germany.

###### Abstract

We present a novel algorithm that allows us to gain detailed insight into the
effects of sparsity in linear and nonlinear optimization. Sparsity is of great
importance in many scientific areas such as image and signal processing,
medical imaging, compressed sensing, and machine learning, as it ensures
robustness against noisy data and yields models that are easier to interpret
due to the small number of relevant terms. It is common practice to enforce
sparsity by adding the $\ell_{1}$-norm as a penalty term. In order to gain a
better understanding and to allow for an informed model selection, we directly
solve the corresponding multiobjective optimization problem (MOP) that arises
when minimizing the main objective and the $\ell_{1}$-norm simultaneously. As
this MOP is in general non-convex for nonlinear objectives, the penalty method
will fail to provide all optimal compromises. To avoid this issue, we present
a continuation method specifically tailored to MOPs with two objective
functions one of which is the $\ell_{1}$-norm. Our method can be seen as a
generalization of homotopy methods for linear regression problems to the
nonlinear case. Several numerical examples – including neural network training
– demonstrate our theoretical findings and the additional insight gained by
this multiobjective approach.

###### Index Terms:

Multiobjective Optimization, Nonsmooth Optimization, Machine Learning,
Sparsity

##  1 Introduction

Optimization problems are at the center of a very large number of problems in
science and engineering. In many situations, the task of minimizing some
specified objective function is further complicated by the desire to enforce
sparsity, i.e., to obtain a solution with as few non-zero entries as possible,
which has several important advantages. Most notably, sparse solutions tend to
be robust against noise. Furthermore, they often allow for a better analysis
and interpretation due to the small number of relevant terms, for instance in
the identification of governing equations [1]. Practical applications can be
found in many different areas, e.g., signal processing, compressed sensing,
neural network training and medical imaging [2, 3, 4, 5, 6].

In signal and image processing, one frequently considers (overdetermined)
linear regression models with an $\ell_{1}$ penalty term:

$\displaystyle\min_{x\in\mathbb{R}^{n}}\left\lVert
Ax-b\right\rVert_{2}^{2}+\lambda\left\lVert x\right\rVert_{1},$ (1)

where $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$ and
$\lambda\in\mathbb{R}^{\geq 0}$. This problem is often referred to as _Lasso_
[7] or _basis pursuit_ [8]. Since appropriately weighting the penalty term is
difficult a priori, the solution is computed for multiple values of $\lambda$,
for instance in _homotopy methods_ [9, 10, 11, 12], where the entire solution
path for $\lambda\in\mathbb{R}^{\geq 0}$ is tracked. The method is based on
the subdifferential of the $\ell_{1}$-norm and the resulting Karush-Kuhn-
Tucker (KKT) conditions, while exploiting the fact that the structure of zero
entries along the solution path is locally constant almost everywhere.

In recent years, nonlinear models have become increasingly popular, in
particular due to the advances in deep learning and the associated neural
network training, where sparsity can help to avoid overfitting [13]. As in the
linear case, the use of $\ell_{1}$ penalty terms is a natural and common
choice for obtaining sparse models, resulting in the problem

$\min_{x\in\mathbb{R}^{n}}f(x)+\lambda\|x\|_{1}\quad\text{with}\quad\lambda\in\mathbb{R}^{\geq
0}.$ (2)

Since $f$ is based on a nonlinear model and is generally non-convex, (2) is
much more difficult to solve than (1).

In contrast to the classical $\ell_{1}$ penalty approach, we will here
interpret sparse optimization as a (non-smooth) _multiobjective optimization
problem_ (MOP), where the objectives $f$ and
$\left\lVert\cdot\right\rVert_{1}$ are minimized simultaneously, i.e.,

$\min_{x\in\mathbb{R}^{n}}\left(\begin{array}[]{c}f(x)\\\ \left\lVert
x\right\rVert_{1}\end{array}\right).$ (MOP-${\ell_{1}}$)

In this context, (2) is equivalent to the _weighted-sum approach_ , i.e.,
$\min_{x\in\mathbb{R}^{n}}\alpha_{1}f(x)+\alpha_{2}\|x\|_{1}$ with
$\alpha_{1},\alpha_{2}\in[0,1]$, $\alpha_{1}+\alpha_{2}=1$ [14]. It is well-
known that every solution of the weighted-sum approach is a solution of the
MOP [14]. As shown in Figure 1, the weight vector $\alpha$ is orthogonal to
the _Pareto front_. Varying $\alpha$ thus yields different compromise
solutions. However, if $f$ is non-convex, then $\alpha$ does no longer
uniquely define a single Pareto optimal point, and we cannot find points in
the non-convex parts of the Pareto front, see [15, 14] for in-depth
discussions. This is also relevant in practice, e.g., in neural network
training. For nonconvex loss functions, well generalizing sparse solutions may
not be computable via the penalty approach. This is demonstrated in Section
4.3.

(a)

(b)

(c)

Figure 1: (a) Different points on the Pareto front by solving
(MOP-${\ell_{1}}$) using the weighted-sum approach with varying $\alpha$. (b)
For non-convex problems, multiple solutions correspond to the same $\alpha$
such that the solution cannot be computed in the non-convex regions. (c) For
strictly concave problems, only the individual minima can be found.

Alternative solution approaches suitable for non-convex problems are
_evolutionary algorithms_ [16, 17], _scalarization methods_ [14], as well as
_continuation methods_ [18, 19, 20]. The latter require sufficiently smooth
objective functions ensuring the existence of a tangent space as the idea is
to compute a discretization of the Pareto critical set by first performing a
so-called _predictor step_ along the tangent space and then computing the
nearest Pareto critical point using an optimization procedure (_corrector
step_).

In this paper, we present a new continuation method for non-smooth MOPs with
two objectives, where the first objective is twice continuously differentiable
and the second objective is the $\ell_{1}$-norm (Section 3). In particular,
the first objective may be non-convex. Although the solution set of non-smooth
MOPs does not generally possess the smooth structure required for classical
continuation, we will show that this specific problem class is piecewise-
smooth. Thus, the solution set can be computed via classical continuation up
to “kinks”. We will show how these kinks can be detected and addressed within
the continuation by exploiting the special structure of the optimality
conditions specifically derived for this purpose. The global minimal point
$0\in\mathbb{R}^{n}$ of the $\ell_{1}$-norm can be used as the starting point.

Compared to scalarization methods (like the _Adam descent method_ [21] applied
to (2)), our approach has the advantage that we compute an approximation of
the entire Pareto set instead of just a single point. In particular, we do not
have to choose a penalty parameter. Compared to evolutionary methods (which
also approximate the entire Pareto set), our approximation is generally closer
to the true Pareto set while requiring fewer function evaluations.
Furthermore, our method detects structural properties like kinks or the
connectedness of the Pareto set, which remain unnoticed in (standard)
evolutionary methods. On the other hand, some of the disadvantages and
challenges of our method will be discussed in Section 5.

Beyond that, alternative approaches at the intersection of sparse and
multiobjective optimization are, for instance, the training of (linear) SVM
models using $\varepsilon$-constraint scalarization [22, 23], or the solution
of MOPs with linear models and the $\ell_{1}$-norm via evolutionary algorithms
[24]. In the context of sparse optimization, generalizations (for convex
problem classes) of the classical homotopy methods were introduced in [25, 26,
27, 28], for instance to arbitrary convex loss-functions or convex penalty
terms. Thus, the main advantage of our method is the treatment of nonlinear,
non-convex problems, which will be demonstrated using several examples of
varying complexity (Section 4).

##  2 Multiobjective and non-smooth optimization

Before presenting the solution approach, we give a short overview of
multiobjective and non-smooth optimization. For more details, see [14, 15, 18]
and [29, 30], respectively.

The general (unconstrained) multiobjective optimization problem is to minimize
multiple (real-valued) functions
$f_{i}:~{}\mathbb{R}^{n}\rightarrow\mathbb{R}$, $i\in\\{1,...,k\\}$, at the
same time:

$\min_{x\in\mathbb{R}^{n}}F(x)=\min_{x\in\mathbb{R}^{n}}\left(\begin{array}[]{c}f_{1}(x)\\\
\vdots\\\ f_{k}(x)\end{array}\right).$ (MOP)

Due to the loss of a total order of $\mathbb{R}^{k}$ for $k>1$, the solution
to (MOP) is not a specific point, but the _set of non-dominated points_ , also
referred to as the _Pareto set_. Non-dominated points $x^{*}$ are
characterized by the fact that there exists no $x$ with $f_{i}(x)\leq
f_{i}(x^{*})$ for all $i\in\\{1,\dots,k\\}$ and $f_{j}(x)<f_{j}(x^{*})$ for at
least one $j\in\\{1,\dots,k\\}$. The image of the Pareto set under $F$ is
known as the _Pareto front_ , lying in the _objective space_ $\mathbb{R}^{k}$.

If all objectives are continuously differentiable, then these points satisfy
the well-known necessary KKT conditions, which state that if a point $x^{*}$
is locally Pareto optimal, there exists a convex combination of the individual
gradients that is zero [18], i.e., there exists
$\alpha^{*}\in(\mathbb{R}^{\geq 0})^{k}$ with

$\displaystyle\sum_{i=1}^{k}\alpha_{i}^{*}\nabla f_{i}(x^{*})=0\text{ and
}\sum_{i=1}^{k}\alpha_{i}^{*}=1.$

Points satisfying the KKT conditions are denoted as _Pareto critical_ , and
the set containing all Pareto critical points is the _Pareto critical set_ ,
which is a superset of the Pareto set.

As the $\ell_{1}$-norm does not satisfy the required smoothness assumptions
for the classical KKT conditions, we need to consider subgradients instead.

######  Definition 2.1 (Subdifferential and subgradient).

Let $f:~{}\mathbb{R}^{n}\rightarrow\mathbb{R}$ be locally Lipschitz
continuous, and let $\Omega\subset\mathbb{R}^{n}$ be the set of points where
$f$ is not differentiable. Then

$\displaystyle\partial f(x):=\mathbf{Conv}(\\{$ $\displaystyle
g\in\mathbb{R}^{n}~{}|~{}\exists(x_{j})_{j}\in\mathbb{R}^{n}\setminus\Omega\text{
with }x_{j}\rightarrow x$ $\displaystyle\text{and }\nabla f(x_{j})\rightarrow
g\text{ for }j\rightarrow\infty\\})$

is the _(Clarke) subdifferential of_ $f$ _in_ $x$, where $\mathbf{Conv}$
denotes the convex hull. An element $g\in\partial f(x)$ is called a
_subgradient_.

If $f$ is continuously differentiable in $x$, then the subdifferential is the
set containing the gradient, i.e., $\partial f(x)=\\{\nabla f(x)\\}$. In case
of the $\ell_{1}$-norm, the subgradients take the following form.

######  Example 2.2 (Subgradients and subdifferential of the $\ell_{1}$-norm).

The $\ell_{1}$-norm in $\mathbb{R}^{n}$, $f_{\ell_{1}}(x)=\left\lVert
x\right\rVert_{1}=\sum_{i=1}^{n}|x_{i}|$, is non-differentiable in all $x$
with $x_{j}=0$ for some $j\in\\{1,\dots,n\\}$. To calculate $\partial
f_{\ell_{1}}(x)$ for an arbitrary $x\in\mathbb{R}^{n}$, we first distinguish
between _active_ and _inactive_ indices, with $\\{j\in J~{}|~{}x_{j}\neq
0\\}=:A(x)$ being the _active set_. It is easy to see that

$g=\left(\begin{array}[]{c}g_{1}\\\ \vdots\\\
g_{n}\end{array}\right)\in\partial
f_{\ell_{1}}(x)\Leftrightarrow\left\\{\begin{array}[]{ll}g_{j}=\textnormal{sgn}(x_{j})&\text{if
}j\in A(x),\\\ g_{j}\in[-1,1]&\text{if }j\notin A(x).\end{array}\right.$

We denote the number of active indices by $n^{A}(x):=\left|A(x)\right|$ and
the number of inactive indices by $n^{0}(x):=n-n^{A}(x)$. The subdifferential
is now the convex hull of $M(x):=2^{n^{0}(x)}$ vectors $g^{i}$ containing the
extreme points of the subdifferential:

$\partial
f_{\ell_{1}}(x)=\mathbf{Conv}\left(\left\\{g^{i}~{}|~{}i\in\\{1,\dots,M(x)\\}\right\\}\right)\text{,
i.e.,}$ $g^{i}_{j}=\textnormal{sgn}(x_{j})~{}\forall j\in
A(x)\quad\text{and}\quad g^{i}_{j}\in\\{-1,1\\}~{}\forall j\notin A(x).$

Consider the point $\hat{x}=(2,0,-3,0)^{\top}\in\mathbb{R}^{4}$ for example.
Then we have $A(\hat{x})=\\{1,3\\}$, $n^{A}(\hat{x})=2$, $n^{0}(\hat{x})=2$,
$M(\hat{x})=4$, and

$\partial f_{\ell_{1}}(\hat{x})=\mathbf{Conv}\left(\left\\{\begin{pmatrix}1\\\
1\\\ -1\\\ -1\end{pmatrix},\begin{pmatrix}1\\\ -1\\\ -1\\\
-1\end{pmatrix},\begin{pmatrix}1\\\ 1\\\ -1\\\
1\end{pmatrix},\begin{pmatrix}1\\\ -1\\\ -1\\\
1\end{pmatrix}\right\\}\right).$

Analogous to the smooth case, necessary optimality conditions for MOPs
containing non-smooth objective functions can be obtained [31].

######  Proposition 2.3 (Non-smooth KKT conditions).

Assume that $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{k}$ is locally Lipschitz
continuous and $x^{*}$ is locally Pareto optimal. Then

$0\in\mathbf{Conv}\left(\bigcup_{i=1}^{k}\partial
f_{i}(x^{*})\right)\subset\mathbb{R}^{n}.$ (3)

Similar to the smooth case, we call $x\in\mathbb{R}^{n}$ _Pareto critical_ if
it satisfies the KKT conditions (3) and we call the set of all those points
the _Pareto critical set_ $\mathcal{P}_{c}$.

##  3 The Continuation Method

Our goal is to compute the entire Pareto critical set $\mathcal{P}_{c}$ of
(MOP-${\ell_{1}}$) using a continuation method that exploits the structure of
the $\ell_{1}$-norm and its subgradients. As a starting point we can use $x=0$
which is always Pareto critical as the minimum of the $\ell_{1}$-norm, and
then move along $\mathcal{P}_{c}$ in the direction where the value of $f$
decreases and $\left\lVert x\right\rVert_{1}$ increases.

Input: $f:\mathbb{R}^{n}\rightarrow\mathbb{R}\in\mathcal{C}^{2}$,
$x_{\text{start}}\in\mathbb{R}^{n}$ satisfying (3)

Output: Discretization of the Pareto critical set
$\mathcal{P}_{\textnormal{PC}}$

$\mathcal{P}_{\textnormal{PC}}=\emptyset$, $x^{0}=x_{\text{start}}$

while End of $\mathcal{P}_{c}$ component has not been reached do

$\mathcal{P}_{\textnormal{PC}}=\mathcal{P}_{\textnormal{PC}}\cup\\{x^{0}\\}$

if $\mathcal{P}_{c}$ is locally smooth then

Predictor step tangential to $\mathcal{P}_{c}$: $x^{p}$ $\leftarrow$
predictor($x^{0}$)

Corrector step with $x^{c}\in\mathcal{P}_{c}$: $x^{c}$ $\leftarrow$
corrector($x^{p}$)

$x^{0}=x^{c}$

else

(De-)Activate entries

end if

end while

Algorithm 1 PC algorithm

Continuation methods (or _predictor-corrector (PC) methods_) yield a sequence
of Pareto critical points by – starting from a Pareto critical point –
performing a prediction step along the tangent space of $\mathcal{P}_{c}$ and
then a correction step which results in the next Pareto critical point close
by, see [18] for details. For the non-smooth problem (MOP-${\ell_{1}}$), the
Pareto critical set may contain kinks such that we cannot exclusively rely on
the tangent space in the prediction step. Thus, when we reach such a kink, we
additionally require a mechanism that computes all possible directions in
which the Pareto critical set might continue, see Algorithm 1 for a sketch. To
this end, we will first have a more detailed look at the KKT conditions (3)
for the special case of problem (MOP-${\ell_{1}}$) in Section 3.1 before
giving a detailed description of the predictor and corrector steps in Sections
3.2 and 3.3. We then discuss strategies for the (de-)activation of indices
(i.e., how to overcome kinks) in Section 3.4, and subsequently address
implementation and globalization in Section 3.5.

###  3.1 Optimality conditions for (MOP-${\ell_{1}}$)

Since we assume $f$ to be at least twice continuously differentiable, the
gradient $\nabla f(x)$ exists for every $x\in\mathbb{R}^{n}$. The
subdifferential of the $\ell_{1}$-norm is given by the convex hull of
$M(x)=2^{n^{0}(x)}$ vectors $g^{i}$ given in Example 2.2. To simplify the
notation, we will write $M$ and $n^{0}$ instead of $M(x)$ and $n^{0}(x)$.
Inserting this into the KKT conditions (3), a point $x$ is Pareto critical for
(MOP-${\ell_{1}}$) if there exist KKT multipliers
$\alpha_{1},\alpha_{2}\in[0,1]$ with $\alpha_{1}+\alpha_{2}=1$ and

$\alpha_{1}\nabla f(x)+\alpha_{2}(\beta_{1}g^{1}+\dots+\beta_{M}g^{M})=0,$ (4)

where $\beta_{i}\in[0,1]$ for $i=1,\dots,M$ and $\sum_{i=1}^{M}\beta_{i}=1$.
The following theorem will further simplify this condition.

######  Theorem 3.1.

Let $x\in\mathbb{R}^{n}$ with $x\neq 0$ and let $a\in
A(x)=\left\\{j\in\\{1,\dots,n\\}~{}|~{}x_{j}\neq 0\right\\}$ be any active
index. Then, the following statements are equivalent:

  1. 1.

$x$ is Pareto critical for (MOP-${\ell_{1}}$)

  2. 2.

The following conditions hold for every $j\in\\{1,\dots,n\\}$:

    1. (a)

If $j\in A(x)$:

      1. i.

$\textnormal{sgn}((\nabla f(x))_{j})=-\textnormal{sgn}(x_{j})$ if $(\nabla
f(x))_{j}\neq 0$ and

      2. ii.

$\left|(\nabla f(x))_{j}\right|=\left|(\nabla f(x))_{a}\right|$

    2. (b)

If $j\notin A(x)$: $\left|(\nabla f(x))_{j}\right|\leq\left|(\nabla
f(x))_{a}\right|$

###### Proof.

We start by showing 1 $\Rightarrow$ 2. Since $x$ is Pareto critical, there
exist $\alpha_{1},\alpha_{2}\in[0,1]$ with $\alpha_{1}+\alpha_{2}=1$ and
$\beta_{i}\in[0,1]$, $i=1,\dots,M$, with $\sum_{i=1}^{M}\beta_{i}=1$ such that
(4) is satisfied.

Note that $\alpha_{1}=0$ would imply that zero is contained in the
subdifferential of the $\ell_{1}$-norm at $x$, which is only the case for
$x=0$. Thus, we must have $\alpha_{1}>0$.

First, assume $j\in A(x)$. Then, $(g^{i}(x))_{j}=\textnormal{sgn}(x_{j})$ for
all $i=1,\dots,M$ and

$\displaystyle 0$ $\displaystyle=\alpha_{1}(\nabla
f(x))_{j}+\alpha_{2}\textnormal{sgn}(x_{j})$ $\displaystyle=\alpha_{1}(\nabla
f(x))_{j}+(1-\alpha_{1})\textnormal{sgn}(x_{j}).$

Since $\alpha_{1}>0$, this is equivalent to

$\displaystyle(\nabla
f(x))_{j}=-\frac{(1-\alpha_{1})}{\alpha_{1}}\textnormal{sgn}(x_{j}).$

Therefore, we have

$\displaystyle\textnormal{sgn}(\nabla
f(x))_{j}=-\textnormal{sgn}(x_{j})\quad\Rightarrow\quad\left|(\nabla
f(x))_{j}\right|=\frac{(1-\alpha_{1})}{\alpha_{1}}.$

Since this has to hold for every active $j$, it follows

$\displaystyle\left|(\nabla f(x))_{j}\right|=\left|(\nabla f(x))_{a}\right|.$

Now assume that $j$ is inactive. Define

$\displaystyle\beta_{j,-}:=\sum_{i:g_{j}^{i}=-1}\beta_{i},\quad\beta_{j,+}:=\sum_{i:g_{j}^{i}=1}\beta_{i}.$

From (4) we obtain

$\displaystyle 0\quad$ $\displaystyle=\alpha_{1}(\nabla
f(x))_{j}+\alpha_{2}(\beta_{j,+}-\beta_{j,-})$
$\displaystyle\Rightarrow\quad\left|(\nabla f(x))_{j}\right|$
$\displaystyle=\left|\beta_{j,-}-\beta_{j,+}\right|\frac{1-\alpha_{1}}{\alpha_{1}}$
$\displaystyle\leq\frac{1-\alpha_{1}}{\alpha_{1}}=\left|(\nabla
f(x))_{a}\right|.$

For the opposite direction ”$\Leftarrow$”, assume that $x$ satisfies 2(a) and
(b). Note that 2(a) implies $\textnormal{sgn}(x_{a})(\nabla f(x))_{a}<0$ and
thus,

$\displaystyle 1-\textnormal{sgn}(x_{a})(\nabla f(x))_{a}>1.$

Define

$\displaystyle\alpha_{1}:=\frac{1}{1-\textnormal{sgn}(x_{a})(\nabla
f(x))_{a}}\in(0,1).$

By 2(a)(i) we obtain

$\displaystyle\alpha_{1}=\frac{1}{1+\left|(\nabla f(x))_{a}\right|}.$

With 2(a)(i) it follows that for all active $j$

$\displaystyle\frac{\alpha_{1}}{1-\alpha_{1}}\left|(\nabla f(x))_{j}\right|$
$\displaystyle=\frac{1}{\frac{1}{\alpha_{1}}-1}\left|(\nabla
f(x))_{a}\right|=1.$

In particular,

$\displaystyle 1=\frac{\alpha_{1}}{1-\alpha_{1}}\left|(\nabla
f(x))_{j}\right|=\frac{\alpha_{1}}{1-\alpha_{1}}\textnormal{sgn}(\nabla
f(x)_{j})\nabla f(x)_{j}$ $\displaystyle\Rightarrow\quad$
$\displaystyle\frac{\alpha_{1}}{1-\alpha_{1}}\nabla
f(x)_{j}=\textnormal{sgn}(\nabla f(x)_{j})=-\textnormal{sgn}(x_{j}).$

Analogously, for all inactive $j$, we can use 2(b) to obtain

$\displaystyle\frac{\alpha_{1}}{1-\alpha_{1}}\left|(\nabla
f(x))_{j}\right|\leq
1\quad\Rightarrow\quad\frac{\alpha_{1}}{1-\alpha_{1}}(\nabla
f(x))_{j}\in[-1,1].$

As a result, the vector $-\frac{\alpha_{1}}{1-\alpha_{1}}\nabla f(x)$ can be
written as a convex combination of the vectors $g^{i}$ defined in Example 2.2.
Let $\beta_{1},...,\beta_{M}$ be the corresponding coefficients and define
$\alpha_{2}=1-\alpha_{1}$. Then

$\displaystyle-\frac{\alpha_{1}}{1-\alpha_{1}}\nabla
f(x)=\beta_{1}g^{1}+\dots+\beta_{M}g^{M}$ $\displaystyle\Leftrightarrow\quad$
$\displaystyle\alpha_{1}\nabla
f(x)+\alpha_{2}(\beta_{1}g^{1}+\dots+\beta_{M}g^{M})=0,$

showing that (4) holds.

∎

Theorem 3.1 states that the absolute values of all gradient entries belonging
to active indices are identical. At the same time, the absolute values of
gradient entries belonging to inactive indices must always be less than or
equal to the absolute value of active gradient entries. Furthermore, the proof
of the theorem yields an explicit formula for the KKT multipliers of Pareto
critical points.

######  Remark 3.2.

If $x$ is Pareto critical for (MOP-${\ell_{1}}$) and $a$ is an active index,
then KKT multipliers corresponding to $x$ are given by

$\displaystyle\alpha_{1}=\frac{1}{1+\left|(\nabla
f(x))_{a}\right|},\quad\alpha_{2}=\frac{\left|(\nabla
f(x))_{a}\right|}{1+\left|(\nabla f(x))_{a}\right|}.$

If $\nabla f(x)\neq 0$ then the KKT multipliers are unique. In particular,
this indicates that “kinks” in the Pareto front can only occur if $\nabla
f(x)=0$.

As mentioned before, since (MOP-${\ell_{1}}$) is a non-smooth problem, the
Pareto critical set will generally contain points in which the tangent space
does not exist. Theorem 3.1 can be used to characterize those points and to
gain insight into the structure of $\mathcal{P}_{c}$ by discriminating between
equality and strict inequality in condition 2(b).

######  Corollary 3.3.

Let $x^{0}\neq 0$ be a Pareto critical point such that the inequality in 2(b)
of Theorem 3.1 is strict. Let $j_{1}<...<j_{n^{A}(x^{0})}$ be the (sorted)
active indices at $x^{0}$. Define the following two mappings that lift the
active components to the $n$-dimensional point $x$ and vice versa:

$\displaystyle h^{1}:\mathbb{R}^{n^{A}(x^{0})}\rightarrow\mathbb{R}^{n},\quad
h^{1}_{j}(\bar{x})=\begin{cases}0&\text{if }j\notin A(x^{0}),\\\
\bar{x}_{i}&\text{if }j_{i}\in A(x^{0}),\end{cases}$ $\displaystyle
h^{2}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n^{A}(x^{0})},\quad
h^{2}_{i}(x)=x_{j_{i}}\quad\forall i\in\\{1,...,n^{A}(x^{0})\\}.$

Then there is an open set $U\subseteq\mathbb{R}^{n}$ with $x^{0}\in U$ such
that $x\in U$ being Pareto critical for (MOP-${\ell_{1}}$) is equivalent to
$A(x)=A(x^{0})$ and $h^{2}(x)$ being Pareto critical for the MOP

$\displaystyle\min_{\bar{x}\in\mathbb{R}^{n^{A}(x^{0})}}\left(\begin{array}[]{c}f(h^{1}(\bar{x}))\\\
\left\lVert\bar{x}\right\rVert_{1}\end{array}\right).$ (7)

###### Proof.

See Section B in the supplementary material. ∎

The previous corollary implies that if the inequality in 2(b) of Theorem 3.1
is strict for a point $x^{0}\in\mathcal{P}_{c}$, then locally around $x^{0}$,
the Pareto critical set of (MOP-${\ell_{1}}$) coincides with the Pareto
critical set of the MOP (7), embedded from $\mathbb{R}^{n^{A}(x^{0})}$ to
$\mathbb{R}^{n}$ via $h^{1}$. The objectives of (7) are smooth around
$h^{2}(x^{0})$ (since all entries of $h^{2}(x^{0})$ are non-zero by
construction), such that we can apply the theory for smooth MOPs to see that
the tangent space of the Pareto critical set $\mathcal{P}_{c}$ must exist in
$x^{0}$. Thus, kinks in $\mathcal{P}_{c}$ can only arise in points where
equality holds in 2(b) for some $j\notin A(x^{0})$. These are precisely the
points where the active set can potentially change, i.e., where an inactive
index might be activated (or vice versa).

Corollary 3.3 shows that it is sufficient to consider the MOP (7) if we are in
a point where the inequality in 2(b) of Theorem 3.1 is strict. Based on (7),
this allows us to construct the predictor (Section 3.2) and corrector (Section
3.3) as in the smooth case except for kinks.

###  3.2 Predictor

In the predictor step, we perform a step along the tangent space of the Pareto
critical set from the starting point $x^{0}\in\mathcal{P}_{c}$. As discussed
above, we assume that an open set around $x^{0}$ exists in which the zero
structure of Pareto critical points does not change, such that the original
MOP locally reduces to the problem (7), where both objective functions are at
least twice continuously differentiable since the $\ell_{1}$-norm is applied
only to non-zero entries. This means that we can compute the predictor
according to classical continuation methods, cf. [18]. More precisely, let
$H:\mathbb{R}^{n^{A}(x^{0})+2}\rightarrow\mathbb{R}^{n^{A}(x^{0})+1}$,

$\displaystyle H(\bar{x},\alpha)=\begin{pmatrix}\alpha_{1}\nabla(f\circ
h^{1})(\bar{x})+\alpha_{2}\nabla f_{\ell_{1}}(\bar{x})\\\
\alpha_{1}+\alpha_{2}-1\end{pmatrix}$

be the map containing the KKT conditions of (7). Then the tangent space of the
Pareto critical set is given by the first $n^{A}(x^{0})$ components of

$\displaystyle\ker(H^{\prime}(h^{2}(x^{0}),\alpha))$ $\displaystyle=$
$\displaystyle\ker\Big{(}\begin{array}[]{c}\alpha_{1}\nabla^{2}(f\circ
h^{1})(h^{2}(x^{0}))\\\ 0\\\ \end{array}$ (10)
$\displaystyle\qquad\qquad\qquad\qquad\begin{array}[]{cc}\nabla(f\circ
h^{1})(h^{2}(x^{0}))&\nabla f_{\ell_{1}}(h^{2}(x^{0}))\\\ 1&1\\\
\end{array}\Big{)}$ (13) $\displaystyle=$
$\displaystyle\ker\left(\begin{matrix}\alpha_{1}\nabla^{2}f(x^{0})_{|A(x^{0})}&\nabla
f(x^{0})_{|A(x^{0})}&s\\\ 0&1&1\\\ \end{matrix}\right),$ (14)

where $(\nabla^{2}f(x))_{|A(x^{0})}$ and $\nabla f(x)_{|A(x^{0})}$ denote the
Hessian and the gradient of $f$ reduced to the active indices, respectively,
and $s=\textnormal{sgn}(h^{2}(x^{0}))=\textnormal{sgn}(x^{0})_{|A(x^{0})}$ is
the result of the componentwise sign function in the active indices of
$x^{0}$. Note that Theorem 3.1 implies $\nabla
f(x^{0})_{|A(x^{0})}=-\left|(\nabla f(x^{0}))_{a}\right|s$.

If the reduced Hessian $\nabla^{2}f(x)_{|A(x^{0})}$ is regular, then the
kernel (i.e., the tangent space) is one dimensional and we obtain a unique
vector (up to scaling and sign). As discussed in [32], we expect this to be
the generic case and assume that this holds for the remainder of this paper.
Otherwise, the kernel could be of higher dimension which would lead to
infinitely many possible directions. In this case, we could try to find a
suitable direction in the kernel by using the directional derivative of the
gradient and the fact that the absolute values of the entries of the gradient
in the active indices must coincide (cf. 2(a) of Theorem 3.1).

If the reduced Hessian $\nabla^{2}f(x)_{|A(x^{0})}$ is regular, then a kernel
vector $\bar{v}$ of (3.2) is given by

$\displaystyle\bar{v}$
$\displaystyle=\gamma\cdot(\bar{v}_{1},\bar{v}_{2},\bar{v}_{3})^{\top}\quad\text{with}\quad\gamma\in\mathbb{R}\quad\text{and}$
(P) $\displaystyle\bar{v}_{1}$ $\displaystyle=\frac{(1+\left|(\nabla
f(x^{0}))_{a}\right|)}{\alpha_{1}}\left(\nabla^{2}f(x^{0})_{|A(x^{0})}\right)^{-1}s,$
$\displaystyle=(1+\left|(\nabla
f(x^{0}))_{a}\right|)^{2}\left(\nabla^{2}f(x^{0})_{|A(x^{0})}\right)^{-1}s,$
$\displaystyle\bar{v}_{2}$ $\displaystyle=1,\qquad\bar{v}_{3}=-1.$

The predictor of (7) is given by $\bar{v}_{1}$, while $\bar{v}_{2}$ and
$\bar{v}_{3}$ represent the corresponding direction in the space of KKT
multipliers. The predictor of (MOP-${\ell_{1}}$) is then
$v_{1}=h^{1}(\bar{v}_{1})$.

####  3.2.1 Direction

If the reduced Hessian is regular, there are only two possible choices for the
direction of the predictor step, namely $v_{1}$ and $-v_{1}$. As we want to
move along the Pareto front in the direction where $f$ becomes smaller (and
the $\ell_{1}$-norm larger), it would be intuitive to choose the sign of
$\nabla f(x)^{\top}v_{1}$ as an indicator, as proposed in [18] and [19]. While
this mostly works, the need to cross a “turning point” may arise. These are
points where we follow the Pareto critical set in the same direction but the
corresponding direction in the objective space vanishes and then turns around.
This phenomenon is not related to the non-smoothness of the $\ell_{1}$-norm
but can also occur in the smooth parts, i.e., on a single _activation
structure_ (the set of points $x\in\mathcal{P}_{c}$ with a constant active set
$A(x)$), as shown in the following example.

######  Example 3.4.

Consider (MOP-${\ell_{1}}$) with

$\displaystyle
f(x)=(x_{1}+1)^{2}+(x_{2}-1)^{4}-\tfrac{1}{2}\left(x_{2}-\tfrac{1}{4}\right)^{3}.$
(15)

Figure 2 shows the Pareto critical set for this problem, which can be computed
by hand in this case. We see that we have two turning points $x^{1}$ and
$x^{2}$ in the Pareto critical set, where the direction in the objective space
turns around.

(a)

(b)

(c)

(d)

Figure 2: (a) $\mathcal{P}_{c}$ for problem (15). (b) Part of
$\mathcal{P}_{c}$ where the turning points $x^{1}$ and $x^{2}$ occur. (c)
Image of $\mathcal{P}_{c}$ with the same coloring as in (b). (d) Derivative of
$f$ in the direction of the predictor, i.e., $\nabla f(x)^{\top}v_{1}$.

We thus use a criterion which stems from classical continuation for dynamical
systems in [32]. It is based on the sign of the determinant of the augmented
jacobian
$\begin{pmatrix}H^{\prime}(h^{2}(x^{0}),\alpha)^{\top},\bar{v}\end{pmatrix}^{\top}$,
which is here given by

$\displaystyle
J_{a}:=\left(\begin{matrix}\alpha_{1}\nabla^{2}f(x^{0})_{|A(x^{0})}&\nabla
f(x^{0})_{|A(x^{0})}&s\\\ 0&1&1\\\
\bar{v}_{1}^{\top}&\bar{v}_{2}&\bar{v}_{3}\\\ \end{matrix}\right).$ (16)

If the sign of $\det{(J_{a})}$ is constant in one activation structure, then
this ensures that $\bar{v}$ has the same orientation with respect to the image
of $H^{\prime}(h^{2}(x^{0}),\alpha)$. Furthermore, it prevents us from walking
back in the direction we came from (cf. [32]). As we will see in Section 3.4,
the correct direction is known at the beginning of a new activation structure.
Therefore, this is a suitable criterion to choose the direction in each
predictor step. Furthermore, a simple calculation shows that

$\displaystyle\det{(J_{a})}$ $\displaystyle=-\bar{v}_{2}\cdot
z\cdot\det{(\alpha_{1}\nabla^{2}f(x^{0})_{|A(x^{0})})},$ (17)

where $z:=2+(1+\left|(\nabla
f(x^{0}))_{a}\right|)^{4}\left\lVert(\nabla^{2}f(x^{0})_{|A(x^{0})})^{-1}s\right\rVert_{2}^{2}>0$
with $a\in A(x^{0})$. As a result, only the determinant of
$\nabla^{2}f(x^{0})_{|A(x^{0})}$ has to be computed in order to determine the
sign of $\det{(J_{a})}$.

####  3.2.2 Step size

Now it remains to choose a step size $h\in\mathbb{R}^{>0}$ that determines how
far we move along the tangent direction, i.e., that determines the predicted
point $x^{p}=x^{0}+h\cdot(\pm v_{1})$. To obtain an even coverage of the
Pareto critical set, we choose a constant step size in the parameter space,
i.e.,

$\displaystyle h:=\frac{\tau}{\left\lVert v_{1}\right\rVert_{2}}$ (18)

with a constant $\tau\in\mathbb{R}^{>0}$. For an even distribution in the
objective space, see [18, 19].

Since we still assume that we are in a part of the Pareto critical set where
the activation structure and sign of the variables does not change, the issue
may occur that $h$ is too large if we choose it according to one of the above
formulas. Therefore, we will choose $\hat{h}=\min(h,\bar{h})$, where $\bar{h}$
determines the maximal length of $v_{1}$, such that the activation structure
and sign does not change, i.e., such that no active index gets inactive or
crosses $0$.

###  3.3 Corrector

The corrector step maps the predicted point onto the Pareto critical set by
finding the nearest point $x^{c}\in\mathcal{P}_{c}$ with the same zero
structure and signs as the predicted point $x^{p}$, i.e., with
$A(x^{c})=A(x^{p})=A(x^{0})$. The standard approach is to solve the system of
equations given by the KKT conditions (4) with $x^{p}$ as starting point.
Since the subgradient of the $\ell_{1}$-norm consists of
$M(x^{0})=2^{n^{0}(x^{0})}$ vectors, this can result in very large systems.
Fortunately, we can exploit Theorem 3.1 to obtain a simpler way, and solve the
following optimization problem:

$\displaystyle\min_{x\in\mathbb{R}^{n}}$ $\displaystyle\left\lVert
x-x^{p}\right\rVert_{2}^{2},\quad s.t.$ (C) $\displaystyle(\nabla
f(x))_{j}^{2}-(\nabla f(x))_{a}^{2}=0$ $\displaystyle\text{ if }x^{0}_{j}\neq
0$ $\displaystyle(\nabla f(x))_{j}^{2}-(\nabla f(x))_{a}^{2}\leq 0$
$\displaystyle\text{ if }x^{0}_{j}=0$ $\displaystyle(\nabla f(x))_{j}\leq 0$
$\displaystyle\text{ if }x^{0}_{j}>0$ $\displaystyle(\nabla f(x))_{j}\geq 0$
$\displaystyle\text{ if }x^{0}_{j}<0$ $\displaystyle x_{j}\geq 0$
$\displaystyle\text{ if }x^{0}_{j}>0$ $\displaystyle x_{j}\leq 0$
$\displaystyle\text{ if }x^{0}_{j}<0$ $\displaystyle x_{j}=0$
$\displaystyle\text{ if }x^{0}_{j}=0$

where $a\in A(x^{0})$ is some index that is currently active.

In practice, we use an SQP method [33] for the solution of (C) with $x^{p}$ as
the initial point. In the case of convergence issues, one can choose a
starting point closer to the Pareto critical point $x^{0}$, i.e.,
$x^{p}+\lambda(x^{0}-x^{p})$ for some $\lambda\in(0,1)$.

###  3.4 Changing the activation structure

So far, we have seen how the predictor and corrector can be computed when the
zero structure is locally constant around $x^{0}$. As discussed at the
beginning of this section, this allows us to compute the smooth parts of the
Pareto critical set in between the kinks. To compute the entire Pareto
critical set (or a connected component of it), a mechanism that is able to
change the activation structure $A(x)$ is required.

According to Theorem 3.1, an index $j$ can only change from inactive to active
(or vice versa) if $x_{j}=0$ and $\left|\nabla f(x)_{j}\right|=\left|\nabla
f(x)_{a}\right|$ with $a\in A(x)$. We will call such an index $j$ _potentially
active_ and denote the set of all potentially active indices by

$\displaystyle A_{p}(x)=\\{j\in\\{1,\dots,n\\}~{}|{}$ $\displaystyle
x_{j}=0\mbox{~{}and}$ (19) $\displaystyle\left|\nabla
f(x)_{j}\right|=\left|\nabla f(x)_{a}\right|\\}.$

####  3.4.1 Activate indices

Let $x^{*}$ be a Pareto critical point with $A_{p}(x^{*})\neq\emptyset$. If
$x^{*}$ is an end point of the current activation structure, then the
construction of the corrector ensures that we reach it after a finite number
of steps (if the step length $\tau$ is sufficiently small). If $x^{*}$ is not
an end point, then we will generally not be able to detect it during
continuation, as we move with a constant step size. Fortunately, the latter
points are less likely to exist, as will be discussed briefly at the end of
this section.

Note that not all potentially active indices necessarily need to be activated,
such that a way to identify the relevant indices is required. Denoting the
number of potentially active indices with $p$, we have $2^{p}$ possible
activation structures. For every possible activation structure
$A^{l}_{p}\subset A_{p}(x^{*})$, $l=1,\dots,2^{p}$, we can compute a vector
$\bar{v}^{l}=(\bar{v}_{1}^{l},\bar{v}_{2}^{l},\bar{v}_{3}^{l})$ as in the
predictor step (Section 3.2). If we assume $\left|\nabla
f(x^{*})_{a}\right|>0$, then according to 2(a) in Theorem 3.1, the correct
sign of the possibly active indices is given by $-\textnormal{sgn}(\nabla
f(x^{*}))_{j}$. Therefore, $\bar{v}^{l}$ is given by

$\ker\left(\begin{matrix}\alpha_{1}\nabla^{2}f(x^{*})_{|A(x^{*})\cup
A^{l}_{p}}&\nabla f(x^{*})_{|A(x^{*})\cup A^{l}_{p}}&s^{l}\\\ 0&1&1\\\
\end{matrix}\right),$ (20)

with $s^{l}=-\textnormal{sgn}(\nabla f(x^{*}))_{|A(x^{*})\cup A^{l}_{p}}$. If
the Hessian $\nabla^{2}f(x^{*})_{|A(x^{*})\cup A^{l}_{p}}$ is regular, we
obtain

$\displaystyle\bar{v}_{1}^{l}$ $\displaystyle=(1+\left|(\nabla
f(x^{*}))_{a}\right|)^{2}\left(\nabla^{2}f(x^{*})_{|A(x^{*})\cup
A^{l}_{p}}\right)^{-1}s^{l},$ (21) $\displaystyle\bar{v}_{2}^{l}$
$\displaystyle=1,$ $\displaystyle\bar{v}_{3}^{l}$ $\displaystyle=-1.$

Analogously to Section 3.2, we obtain the direction
$v_{1}^{l}\in\mathbb{R}^{n}$ by the corresponding embedding
$h^{1}:\mathbb{R}^{n^{A}(x^{*})+\left|A_{p}^{l}\right|}\rightarrow\mathbb{R}^{n}$.
If $A^{l}_{p}$ is a set of potentially active indices that can actually be
activated, i.e., if there is a curve $\gamma$ of Pareto critical points
starting in $x^{*}$ with the active set $A(x^{*})\cup A^{l}_{p}$, then by
construction, $v_{1}^{l}$ is the tangent vector of that curve. By analyzing
the relationship between the derivative of $\gamma$ at $x^{*}$ and the
conditions in Theorem 3.1, we obtain necessary conditions for $A^{l}_{p}$ to
be activated. This is done in the following lemma and the subsequent
corollary.

######  Lemma 3.5.

Let $x^{*}\in P_{c}\setminus\\{0\\}$ with $\nabla f(x^{*})\neq 0$. Let
$\gamma:(-1,1)\rightarrow\mathbb{R}^{n}$ such that

  * (i)

$\gamma$ is continuously differentiable,

  * (ii)

$\gamma(t)\in P_{c}\quad\forall t\in[0,1)$,

  * (iii)

$\gamma(0)=x^{*}$,

  * (iv)

$A(\gamma(t))=:A_{\gamma}$ is constant $\forall t\in(0,1)$,

Let $a\in A(x^{*})$ and $w=\gamma^{\prime}(0)$. Then for all $j\in
A_{\gamma}\setminus A(x^{*})$ we have

  1. (a)

either $w_{j}=0$ or $\textnormal{sgn}(w_{j})=-\textnormal{sgn}((\nabla
f(x^{*}))_{j})$,

and for all $j\in A_{p}(x^{*})\setminus A_{\gamma}$ we get

  1. (b)

$\textnormal{sgn}((\nabla f(x^{*}))_{j})(\nabla^{2}f(x^{*}))_{j}w$
$\leq\textnormal{sgn}((\nabla f(x^{*}))_{a})(\nabla^{2}f(x^{*}))_{a}w$,

where $(\nabla^{2}f(x^{*}))_{j}$ is the $j$-th row of $\nabla^{2}f(x^{*})$.

###### Proof.

See Section C in the supplementary material. ∎

Unfortunately, without knowledge of the Pareto critical set, the existence of
the curve $\gamma$ in the previous lemma is almost impossible to verify in
practice. Recall that Corollary 3.3 indicates that $\mathcal{P}_{c}$ is
piecewise smooth, where parts in between kinks are embeddings of Pareto
critical sets of the smooth problems (7) into $\mathbb{R}^{n}$. Thus, the
existence of $\gamma_{|(0,1)}$ follows from the smooth structure of (7). In
light of this, the existence of $\gamma$ implies a regularity of the smooth
parts when they reach a kink. For example, roughly speaking, continuous
differentiability of $\gamma$ on an open neighborhood of $t=0$ implies that
the “curvature” of $\mathcal{P}_{c}$ is bounded.

Note that the statement of the previous lemma only involves the derivative
$w=\gamma^{\prime}(0)$ and the new active set $A_{\gamma}$. By the
construction in (21), $w$ is (up to scaling) equal to the direction
$v_{1}^{l}$ for the index set $A^{l}_{p}=A_{\gamma}\setminus A(x^{*})$.
Although $A_{\gamma}$ is not known in practice, we can use this relation to
identify the indices that need to be activated, as shown in the following
corollary.

######  Corollary 3.6.

Let $x^{*}\in P_{c}\setminus\\{0\\}$ with $\nabla f(x^{*})\neq 0$, let
$\gamma$ be a curve as in Lemma 3.5 and let $a\in A(x^{*})$.
Let $v^{l}_{1}$ be the vector computed via (21) and the appropriate embedding
$h^{1}$ corresponding to the index set $A^{l}_{p}=A_{\gamma}\setminus
A(x^{*})$. Then, for $w=\sigma v_{1}^{l}$ with $\sigma\in\\{-1,+1\\}$,

  1. (a)

either $w_{j}=0$ or $\textnormal{sgn}(w_{j})=-\textnormal{sgn}((\nabla
f(x^{*}))_{j})$,

and for all $j\in A_{p}(x^{*})\setminus A_{\gamma}$,

  1. (b)

$\textnormal{sgn}((\nabla
f(x^{*}))_{j})(\nabla^{2}f(x^{*}))_{j}w\leq-\sigma(1+\left|(\nabla
f(x))_{a}\right|)^{2}$,

where $(\nabla^{2}f(x^{*}))_{j}$ is the $j$-th row of $\nabla^{2}f(x^{*})$.

######  Remark 3.7.

According to (a) in Corollary 3.6, for all possible activation structures
$A_{p}^{l}(x^{*})\neq\emptyset$, there either exists only one valid direction
($+v_{1}^{l}$ or $-v^{l}_{1}$) or the tangent vector $w$ of the corresponding
curve $\gamma$ has to satisfy $w_{j}=0$ for all $j\in A_{p}(x^{*})$, i.e.,
$\gamma$ has to be tangential to a vector with activation structure
$A(x^{*})$, which is unlikely.

Using the conditions from Corollary 3.6, we now obtain a predictor for the
suitable activation structures with the correct sign. Since we can exclude the
direction we are coming from, we have to follow each of the remaining
directions. This is shown in Examples 3.8 and 3.9.

######  Example 3.8.

Solve (MOP-${\ell_{1}}$) with

$\displaystyle f(x)=(x_{1}-2)^{2}+(x_{2}-1)^{2}+(x_{3}-1)^{2}.$ (22)

Figure 3(a) shows the Pareto critical set for this problem. Consider
$x^{*}=(1,0,0)^{\top}$. Since $\nabla f(x^{*})=(-2,-2,-2)^{\top}$, we have
$A_{p}(x^{*})=\\{2,3\\}$ (cf. (19)). According to (21), we get the following
possible directions:

$\displaystyle v_{1}^{0}=\begin{pmatrix}\frac{9}{2}\\\\[3.00003pt]
0\\\\[3.00003pt]
0\end{pmatrix},v_{1}^{1}=\begin{pmatrix}\frac{9}{2}\\\\[3.00003pt]
\frac{9}{2}\\\\[3.00003pt]
0\end{pmatrix},v_{1}^{2}=\begin{pmatrix}\frac{9}{2}\\\\[3.00003pt]
0\\\\[3.00003pt]
\frac{9}{2}\end{pmatrix},v_{1}^{3}=\begin{pmatrix}\frac{9}{4}\\\\[3.00003pt]
\frac{9}{4}\\\\[3.00003pt] \frac{9}{4}\end{pmatrix}$

After checking condition (a) in Corollary 3.6, we know that only

$\displaystyle v_{1}^{0},-v_{1}^{0},v_{1}^{1},v_{1}^{2}\text{ and }v_{1}^{3}$

remain as suitable directions. By checking condition (b), the set of
admissible directions can be reduced to

$\displaystyle-v_{1}^{0}\text{ and }v_{1}^{3}.$

Assuming that we started in the origin, we have to walk in direction
$v_{1}^{3}$, i.e., both $j=2$ and $j=3$ become active (and positive).

######  Example 3.9.

Solve (MOP-${\ell_{1}}$) with

$\displaystyle f(x)=(x_{1}-2)^{2}+(x_{2}-1)^{2}+2x_{1}x_{3}.$ (23)

Figure 3(b) shows the relevant part of $\mathcal{P}_{c}$ for this problem.
Consider $x^{*}=(1,0,0)^{\top}$. Since $\nabla f(x^{*})=(-2,-2,2)^{\top}$, we
have $A_{p}(x^{*})=\\{2,3\\}$. According to (21), we get as possible
directions:

$\displaystyle v_{1}^{0}=\begin{pmatrix}\frac{9}{2}\\\\[3.00003pt]
0\\\\[3.00003pt]
0\end{pmatrix},v_{1}^{1}=\begin{pmatrix}\frac{9}{2}\\\\[3.00003pt]
\frac{9}{2}\\\\[3.00003pt]
0\end{pmatrix},v_{1}^{2}=\begin{pmatrix}-\frac{9}{2}\\\\[3.00003pt]
0\\\\[3.00003pt]
9\end{pmatrix},v_{1}^{3}=\begin{pmatrix}-\frac{9}{2}\\\\[3.00003pt]
\frac{9}{2}\\\\[3.00003pt] 9\end{pmatrix}$

Verification of conditions (a) and (b) in Corollary 3.6 reduces the set of
admissible directions to $-v_{1}^{0}$, $v_{1}^{1}$ and $v_{1}^{2}$. Assuming
that we started in the origin, we have to (try to) advance with the
continuation into the directions $v_{1}^{1}$ and $-v_{1}^{2}$. A predictor
step along $v_{1}^{1}$ directly terminates since the corrector delivers
$x^{*}$ as the nearest Pareto critical point and therefore, the remaining
direction is $-v_{1}^{2}$. Thus, only $j=3$ is activated.

(a)

(b)

Figure 3: (a) $\mathcal{P}_{c}$ for problem (22). (b) Part of
$\mathcal{P}_{c}$ of problem (23).

In practice, often only two possible directions remain after applying
Corollary 3.6 – the old and new direction. Nevertheless, it can occur that the
Pareto set splits, i.e., that there are at least 3 suitable directions. In
this case, the continuation method needs to be advanced in each direction
individually. An example where this is demonstrated can be found in Section A
of the supplementary material. If no direction satisfies the conditions in
Corollary 3.6, we have reached the end of a connected component of
$\mathcal{P}_{c}$. In the smooth parts of the Pareto critical set with
constant activation structure, a connected component can only end if $\nabla
f(x)=0$, i.e., a local minimizer for the main objective has been found [34].

######  Remark 3.10.

Although the requirements of Corollary 3.6 do not hold if $\nabla f(x)=0$,
similar conditions can be obtained for this case. First, note that $\nabla
f(x)=0$ implies that $A_{p}(x^{*})$ is the set of all inactive indices. Since
we cannot predict the sign of $x_{j}$ with $j\in A_{p}(x^{*})$ via the sign of
the gradient, we have to consider every possible combination of the signs of
potentially active indices for every possible activation structure
$A_{p}^{l}$. The conditions of Corollary 3.6 can then be used as necessary
conditions, where the sign of the gradient is replaced by the opposite sign of
the currently considered activation and sign structure. An additional
condition can be obtained by using the fact that $\textnormal{sgn}(\nabla
f(x)_{j})=-\textnormal{sgn}(x_{j})$ for all $j\in A(x)$ has to hold. If the
number of inactive indices is large, this results in a very large number of
directions that have to be checked. Therefore a more efficient approach needs
to be developed in future work for very high-dimensional problems.

Finally, we note that it is in theory possible that indices are activated in
points besides the end points of activation structures. In other words, the
Pareto critical set may split into two branches, with one branch continuing
along the old structure. Since we move through smooth parts of the Pareto
critical set with discrete steps, points like this cannot be detected by our
method. But based on Corollary 3.6 and Remark 3.7, one can show that such
points are unlikely to exist in practice.

####  3.4.2 Deactivate indices

A situation where indices have to be deactivated is when the predictor crosses
zero for some entry, i.e.,
$\textnormal{sgn}(x^{p}_{j})\neq\textnormal{sgn}(x^{0}_{j})$ for some $j\in
A(x^{0})$. In this situation, we reduce the predictor step length $h$ to
$\bar{h}$ as the largest step size such that

$\displaystyle x^{p}_{j}\leq 0\ \text{if}\ x^{0}_{j}<0\quad\mbox{and}\quad
x^{p}_{j}\geq 0\ \text{if}\ x^{0}_{j}>0,$

for all $j\in A(x^{0})$, and the indices in $A(x^{0})\setminus A(x^{p})$ are
considered as inactive in the subsequent corrector step. If the corresponding
corrector step (C) has a solution, then the continuation method can continue
on the new activation structure. Otherwise, we choose a step length in
$(0,\bar{h})$ to find a new Pareto critical point with respect to the old
active set $A(x^{0})$. The likelihood of an active index being falsely
deactivated by this mechanism depends on the constant step size $\tau$ and the
curvature of the Pareto critical set.

A second scenario where active indices can be deactivated is the case where
the corrector computes a point $x^{c}$ where one of the previously active
indices is now zero, i.e., $x^{c}_{i}=0$ for some $i\in A(x^{0})$. This
implies that some of the previously active indices $A(x^{0})$ are now only
potentially active. In this case, we proceed as in the activation procedure
described in Section 3.4.1.

###  3.5 Implementation & Globalization

We now have all the ingredients to implement Algorithm 1.
Activation of indices:

  1. 1.

If $A_{p}(x)$ (Eq. (19)) is non-empty, calculate $v^{l}$ for all
$A_{p}^{l}\subset A_{p}(x)$ using (21),

  2. 2.

Reduce number of directions using Corollary 3.6,

  3. 3.

If more than one potential direction remains, proceed in each direction. In
most cases, the corrector will yield the information that $\mathcal{P}_{c}$
only continues in one direction.

Predictor step (when $A_{p}(x)=\emptyset$):

  1. 1.

Compute direction $v_{1}=h^{1}(\bar{v}_{1})$ using (P),

  2. 2.

Determine sign of $v_{1}$ using (17),

  3. 3.

Compute the step size along $\pm v_{1}$ using (18),

  4. 4.

Determine $x^{p}=x^{0}+h\cdot(\pm v_{1})$.

Deactivation of indices: Set the indices $j$ to zero for which
$\textnormal{sgn}(x^{p}_{j})$ changes and reduce the predictor step length $h$
to $\bar{h}$ to avoid sign changes.
Corrector step: Compute point $x^{c}\in\mathcal{P}_{c}$ using (C) with $x^{p}$
as the initial condition.

For the starting point, an intuitive choice would be $x_{\text{start}}=0$,
which is Pareto critical and even Pareto optimal as the global minimal point
of the $\ell_{1}$-norm. According to Theorem 3.1, the first index that has to
be activated in $x_{\text{start}}=0$ is the one which corresponds to the
largest absolute value of the gradient. If there are more than one maximal
entries, we can apply the strategy described in Section 3.4. Of course, every
other Pareto critical point would also work as a starting point, from which we
would have to walk in both directions of the predictor.

By construction, our method only computes the connected component of the
Pareto critical set that contains the starting point $x_{\text{start}}$. Thus,
if there are multiple connected components, then additional starting points
are required to restart our method once the first component is computed. In
practice, this can be done by considering the objective space: Due to the
definition of Pareto optimality, the new starting point should have a smaller
$f$ value and a larger $\ell_{1}$-norm or vice versa (depending on if we
reached the lower end or the upper end of the image of the component). This
defines an area in the objective space in which the image of the new point
should lie. To actually find it, the $\varepsilon$_-constraint method_ [15]
can be used, which is a solution method that can be restricted to certain
parts of the objective space. Alternatively, _reference point methods_ [15]
can be used, where the solution can be influenced by choosing an appropriate
_reference point_ in the objective space.

Finally, it is important to note that our method computes points that are
Pareto _critical_ , but not necessarily Pareto _optimal_. To obtain the actual
Pareto set, we can apply a _non-dominance test_ [35]. Since the output of our
method is a fine pointwise approximation of the Pareto critical set (which
contains the Pareto set), this will generally result in a fine approximation
of the Pareto set.

##  4 Numerical Results

We now study the behavior of the continuation method using three examples. We
begin with two nonlinear toy examples, and then address the identification of
governing equations from data via the SINDy algorithm [1]. Note that the
latter is a regression problem and can consequently be adressed by existing
homotopy approaches as well as the weighted sum method. Finally, the training
of a neural network is studied. While the cost of computing the entire Pareto
set via continuation is certainly higher in comparison to calculating only a
single solution, we note that (a) the individual corrector problems can be
solved rather quickly due to good initial guesses and (b) the randomness (and
thus, the necessity to perform multiple runs) can be reduced.

###  4.1 Toy-Examples

As a first example and to highlight the important aspects of nonlinear
continuation, we consider (MOP-${\ell_{1}}$) with

$f(x)=\left(x_{1}-\tfrac{1}{4}\right)^{2}+\left(x_{2}-\tfrac{1}{2}\right)^{2}+\left(x_{3}-1\right)^{4}-\tfrac{1}{2}\left(x_{3}-\tfrac{1}{4}\right)^{3}.$

Figure 4 shows the result of our method.

(a)

(b)

(c)

(d)

Figure 4: (a) $\mathcal{P}_{c}$ for the toy example. (b) Image of
$\mathcal{P}_{c}$. (c) and (d) Comparison to the weighted sum method, with
which the red parts cannot be computed.

We observe that indices are activated and deactivated multiple times. Points
where the activation structure changes are $x^{1}=(0,0,0.375)$,
$x^{2}=(0,0.25,\frac{83}{151})$, $x^{3}\approx(0,0.25,0.81)$,
$x^{4}\approx(0,0,1.07)$, $x^{5}\approx(0,0,1.93)$ and
$x^{6}\approx(0,0.25,2.014)$. In Figure 4(c) and (d), we see that only the
convex part of $\mathcal{P}_{c}$ can be computed with the weighted sum method.
In particular, a large part where only $j=3$ is active is missing.

To demonstrate the feasibility in high dimensions, let

$f(x)=(x_{1}-a_{1})^{4}+\sum_{i=2}^{n}(x_{i}-a_{i})^{2},$ (24)

where $a_{1}=1$ and $a_{i},i=2,\dots,n$, are random values in $[-1,1]$
(similar to [19], Problem 4). Figure 5 shows the result of our method for
$n=1000$ with step size $\tau=0.25$, which coincides with the analytical
solution.

Figure 5: Pareto front for (MOP-${\ell_{1}}$) with $f$ according to (24).

###  4.2 SINDy

The second example addresses the sparse identification of nonlinear dynamical
systems [1], where the governing system of ordinary (or partial) differential
equations is sought:

$\dot{x}(t)=g(x(t)).$

Here, $x(t)\in\mathbb{R}^{m}$ is the $m$-dimensional system state at time $t$,
whereas the function $g$ describing the system dynamics is unkown. The SINDy
approach is then to represent $g$ via a dictionary of $c$ Ansatz functions
(denoted by $\Theta(x)$) with the corresponding coefficients
$W\in\mathbb{R}^{c\times m}$:

$g(x(t))=W^{\top}\Theta(x(t)).$

The entries of $\Theta$ are arbitrary linear and nonlinear functions such as
sine and cosine or polynomials. Using $N$ (potentially noisy) training data
points $((x^{1},\dot{x}^{1}),...,(x^{N},\dot{x}^{N}))$, the objective is to
identify a sparse representation of $g$ in terms of $\Theta$, and thus to
obtain a robust and interpretable dynamical system from data. The
corresponding main objective is thus the following least squares term:

$f(W)=\frac{1}{mcN}\sum_{i=1}^{N}\left\lVert
W^{\top}\Theta(x^{i})-\dot{x}^{i}\right\rVert_{2}^{2}.$

As this is a convex function, the corresponding MOP is also convex and can
thus be addressed using the standard penalization approach, i.e., the weighted
sum method $\min_{W}f(W)+\lambda\|W\|_{1}$ (where $W$ is transformed to one
large weight vector). However, as we will see, the solution may be highly
sensitive to the choice of $\lambda$.

We consider the well-known Lorenz system with parameters $\sigma=10$,
$\rho=28$ and $\beta=\frac{8}{3}$ for which the dynamical system posses a
chaotic solution:

$\dot{x}=\begin{pmatrix}\sigma\cdot(x_{2}-x_{1})\\\ x_{1}\cdot(\rho-
x_{3})-x_{2}\\\ x_{1}\cdot x_{2}-\beta\cdot x_{3}\end{pmatrix}.$

We collect $N=10,000$ training data points of a single long-term simulation to
which we add white noise. The dictionary $\Theta$ consists of polynomial terms
up to order $3$ and is thus comprised of $c=20$ terms, i.e.,
$W\in\mathbb{R}^{20\times 3}$.

(a)

(b)

(c)

Figure 6: (a) Convex Pareto front for the SINDy approach applied to the Lorenz
system. (b) Evolution of all $W_{i,j}$ which become greater than 0.5 over the
iterations. (c) Trajectory of the Lorenz system (black) and the resulting
trajectories of the approximated solutions where the corresponding point on
the Pareto front is marked in (a).

The Pareto front of the corresponding MOP – obtained via continuation – is
shown in Figure 6(a), the corresponding non-zero weights and three
trajectories are visualized in (b) and (c), respectively. We see that although
the Pareto front is convex, choosing the weight in a penalization approach can
be challenging. Large parts of the Pareto front are almost parallel to either
the horizontal or vertical axis, such that a small change in the penalty
weight will result in large jumps along the front (cf. also Figure 1). Note
that the so-called _knee_ point is not necessarily the desired compromise, as
the system corresponding to the turquoise point is significantly less accurate
on unseen data than the solution corresponding to the yellow point.

###  4.3 Neural Network

Our next example is the training process of the weights $W$ of a feed-forward
neural network (NN), i.e., $y=g_{W}(x)$ [13]:

$\displaystyle z_{0}$ $\displaystyle=x,$ $\displaystyle z_{j+1}$
$\displaystyle=h\left(w_{j}^{\top}z_{j}+w_{j}^{0}\right)\quad\mbox{for
}j=0,\ldots,K-1,$ $\displaystyle y$ $\displaystyle=\sigma(z_{K}),$

where $K$ is the number of _hidden layers_ and $\sigma$ and $h$ are activation
functions. The set of weight matrices, i.e.,
$W=\\{(w_{1}^{0},w_{1})\ldots,(w_{K-1}^{0},w_{K-1})\\}$ (which can again be
transformed into a single weight vector), is chosen such that the error on a
given training data set $((x^{1},y^{1}),\ldots,(x^{N},y^{N}))$ is minimized.
In the case of a classification problem, the error is expressed as the mean of
the cross entropy loss, which results in the following objective function:

$f(W)=-\frac{1}{N}\sum_{i=1}^{N}\sum_{c=1}^{M}(y^{i}_{c}\log(g_{W}(x^{i})_{c})),$

where $M$ is the number of classes in the data set and $y_{c}^{i}\in\\{0,1\\}$
indicates whether $x^{i}$ belongs to class $c$. In addition, sparsity is often
desired to increase robustness against noise.

We here consider a feed forward neural network with $K=2$ hidden layers to
solve the classification task of the well-known Iris data set [36]. The
activation function in the hidden layers is $h(x)=\tanh(x)$, and we chose the
softmax function for the output layer, i.e.,
$\sigma(x)_{i}=\frac{\exp(x)_{i}}{\sum_{j=1}^{m}(\exp(x)_{j})}$. Each layer
consists of two neurons such that the total number of weights is $|W|=25$.

Figure 7: Pareto front for the NN training on the Iris data set (black). The
colored dots show solutions obtained via Adam with different weights for the
$\ell_{1}$ penalty term. These cluster in two places, cf. the zoomed in plots
in the bottom.

Figure 7 shows the Pareto front obtained via the continuation method when
starting in $W=0$ and taking $\tau=0.1$ as step size. We computed three
different connected components of the Pareto critical set and applied a non-
dominance test afterwards. Since the images of the connected components
intersect, they do not appear as distinct connected components in the
objective space. Note that large parts of the front are non-convex. To compare
the solution to stochastic gradient descent training, we also show several
solutions obtained using the Adam descent method [21] with 5000 epochs and
random initial conditions as well as varying weights $\lambda$ for the penalty
term $\|W\|_{1}$. When choosing $\lambda$ too large, the solution always tends
towards zero (which is likely due to the non-convexity and the convergence to
a local optimum). For small values of $\lambda$ (including $\lambda=0$), we
obtain points close to the Pareto front. However, these have larger values of
$\|W\|_{1}$ for a comparable training error. More critically, for
$\lambda=10^{-3}$ we obtain solutions in both clusters, depending on the
initial condition. Finally, solutions with intermediate sparsity are never
obtained using Adam, as these appear to be – in the dynamical systems sense –
less attractive for gradient-descent schemes. This demonstrates the
difficulties of the standard penalization approach and the clear advantage of
the continuation method, which produces high-quality trade-off solutions right
away. This is further emphasized when considering the test error based on data
that was not used for training. Figure 8 shows the test error for the
solutions computed with the continuation method and Adam. While both
approaches are capable of computing solutions with a low training error, we
observe that the larger sparsity of the continuation method is significantly
more robust in terms of the test error.

Figure 8: Test error (black) and training error (grey) against $\ell_{1}$-norm
for Pareto critical points computed via continuation. The colored dots show
the test error of the Adam solutions.

##  5 Towards high-dimensional problems

Except for the polynomial toy example, the experiments presented in the last
section were of low to moderate dimensions. In this section, we discuss the
general scalability of our method towards higher dimensions and discuss
several drawbacks and approaches to resolve some of these. The greatest
challenges are related to symmetries, which is a well-known problem in
continuation methods [32] as well as in other training approaches for neural
networks [37].
Open problems:

  1. 1.

The computation of the (reduced) Hessian and its inverse are very expensive
with increasing dimension $n$.

  2. 2.

The necessary conditions for the activation of indices introduced in Section
3.4.1 often yield a unique direction in which to continue, but there are also
cases where several directions remain and the Pareto critical set splits into
multiple branches, see also the example in Section A of the supplementary
material. In particular, this is the case if $f$ possesses symmetries, as is
often the case for neural networks.

  3. 3.

The Hessian may become singular in particular situations, for instance when
symmetries occur.

  4. 4.

With increasing dimension, the $\varepsilon$-constraint method which we use to
find a new component, becomes more expensive and sometimes numerically
unstable.

In order to address issue 1), we can use a BFGS or SR1 update (cf. [33])
during the optimization in the Corrector step to obtain an approximation of
the Hessian for the predictor step. Furthermore, this way we can directly
compute the inverse by a rank-1-update. Experiments have shown that the SR1
update is more suitable in our case. To handle problem 2), we can apply a
greedy strategy, i.e., we proceed only in the first suitable direction and
ignore the others. This is reasonable in particular for splitting caused by
symmetry, as the branches of $\mathcal{P}_{c}$ then correspond to the same
part of the Pareto front. This way, we can also avoid most of the cases where
the Hessian is singular (issue 3)). Nevertheless, it is an open problem how to
continue when the Hessian matrix is not regular. The only way to solve problem
4) is to not switch the component. Therefore, one could compute only the
_relevant_ connected component by computing the appropriate
$x_{\textnormal{start}}$ using scalarization (e.g., weighted sum solved with
Adam). Alternatively, additional knowledge (e.g., regarding the influence of
bias neurons on the number of connected components in a neural network) may
help in finding suitable transitions to neighboring connected components. In
the following example, we apply the above suggestions to train a larger neural
network.

######  Example 5.1 (NN for reduced MNIST).

As in Section 4.3, we train a feed-forward neural network. We consider the
well-known MNIST data set [38]. The data set consists of images with $28\times
28$ pixels which results in a large number of trainable parameters even for
few neurons. Therefore, we only take images labeled with ”$3$” or ”$6$” and
reduce the number of pixels to $6\times 6$, similar to [39]. To classify these
images, we use a neural network with $K=2$ hidden layers and four neurons per
layer, which results in $\left|W\right|=173$ parameters to train. We choose
the softplus activation function $h(x)=\log(\exp(x)+1)$. Figure 9(a) shows the
Pareto front obtained via the continuation method when starting in $W=0$ with
$\tau=0.25$. Similar to Section 4.3, we computed three components and used a
non-dominance test to obtain the results. The shown Adam solutions are
obtained with $500$ epochs.

(a)

(b)

Figure 9: (a) Pareto front for the NN training on the reduced MNIST data set
(black). The colored dots show solutions obtained via Adam with different
weights for the $\ell_{1}$-term. (b) Test error (black) and training error
(grey) against the $\ell_{1}$-norm for Pareto critical points computed via
continuation. The colored dots show the test error of the Adam solutions.

We observe again that they cluster in two places and that they are
significantly less sparse. Furthermore, the Adam solutions suggest that there
exists another part of the Pareto front with significantly less sparsity but a
slightly lower error, which is missing in our solution and would require
additional runs with additional $x_{\text{start}}$. However, as can be seen in
Figure 9(b), the effect on both training and test error is almost negligible
and thus, these solutions are not of particular practical relevance. Observe
that the test error does not increase again, as the best fitting model is not
sparse.

##  6 Conclusion

The continuation method we present in this paper allows researchers and
practitioners with an interest in sparse optimization to systematically
calculate the entire set of compromise solutions between sparsity and the main
objective, thereby allowing for much more insight as well as informed model
selection than weighted-sum approaches or descent methods providing single
Pareto optima [21, 40]. The presented approach extends the well-known homotopy
methods to the much more complex nonlinear problem setting, where weighting
approaches fail due to non-convexity. The optimality conditions specifically
derived for this problem class allow us to use continuation methods from
smooth multiobjective optimization almost everywhere, and to efficiently
analyze and exploit the problem structure in points where non-smoothness
occurs.

To unfold the full potential of this multiobjective approach to sparse
optimization, several additional questions should be addressed in future work.
These contain, for instance, how to deal with symmetries, in particular with
regard to higher-dimensional problems. This is of particular interest for the
training of neural networks with many training parameters. Furthermore, in
order to further extend the range of possible applications, the consideration
of additional – more general non-smooth – objectives as well as constraints
(see, for instance, [41]) is certainly of high interest. Finally, an extension
to more than one smooth objective may be useful. This could be done by
combining our results with the approach in [20].

## Acknowledgments

This research has been funded by the European Union and the German Federal
State of North Rhine-Westphalia within the EFRE.NRW project “SET CPS”, and by
the DFG Priority Programme 1962 “Non-smooth and Complementarity-based
Distributed Parameter Systems”.

## References

  * [1] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [2] Y. C. Eldar and G. Kutyniok, _Compressed Sensing: Theory and applications_. Cambridge University Press, 2012.
  * [3] D. Vidaurre, C. Bielza, and P. Larrañaga, “A survey of L1 regression,” _International Statistical Review_ , vol. 81, no. 3, pp. 361–387, 2013.
  * [4] R. G. Baraniuk, E. Candes, M. Elad, and Y. Ma, “Applications of sparse representation and compressive sensing,” _Proceedings of the IEEE_ , vol. 98, no. 6, pp. 906–909, 2010.
  * [5] M. Elad, M. A. T. Figueiredo, and Y. Ma, “On the role of sparse and redundant representations in image processing,” _Proceedings of the IEEE_ , vol. 98, no. 6, pp. 972–982, 2010.
  * [6] M. D. Plumbley, T. Blumensath, L. Daudet, R. Gribonval, and M. E. Davies, “Sparse representations in audio and music: From coding to source separation,” _Proceedings of the IEEE_ , vol. 98, no. 6, pp. 995–1005, 2010\.
  * [7] R. Tibshirani, “Regression shrinkage and selection via the Lasso,” _Journal of the Royal Statistical Society. Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996.
  * [8] S. Chen and D. Donoho, “Basis pursuit,” in _Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers_ , vol. 1, 1994, pp. 41–44.
  * [9] M. R. Osborne, B. Presnell, and B. A. Turlach, “A new approach to variable selection in least squares problems,” _IMA Journal of Numerical Analysis_ , vol. 20, no. 3, pp. 389–403, 2000.
  * [10] D. M. Malioutov, M. Cetin, and A. S. Willsky, “Homotopy continuation for sparse signal representation,” in _Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP’05)_ , vol. 5. IEEE, 2005, pp. 733–736.
  * [11] B. Bringmann, D. Cremers, F. Krahmer, and M. Möller, “The homotopy method revisited: Computing solution paths of $\ell_{1}$-regularized problems,” _Mathematics of Computation_ , vol. 87, no. 313, pp. 2343–2364, 2018.
  * [12] D. L. Donoho and Y. Tsaig, “Fast solution of $\ell_{1}$-norm minimization problems when the solution may be sparse,” _IEEE Transactions on Information Theory_ , vol. 54, no. 11, pp. 4789–4812, 2008.
  * [13] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_. MIT Press, 2016, http://www.deeplearningbook.org.
  * [14] K. Miettinen, _Nonlinear multiobjective optimization_ , 4th ed., ser. International series in operations research & management science. Kluwer Academic Publishing, 2004.
  * [15] M. Ehrgott, _Multicriteria optimization_ , 2nd ed. Springer, 2005.
  * [16] K. Deb, _Multi-objective optimization using evolutionary algorithms_. John Wiley & Sons, 2001, vol. 16.
  * [17] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” _IEEE transactions on evolutionary computation_ , vol. 6, no. 2, pp. 182–197, 2002.
  * [18] C. Hillermeier, _Nonlinear multiobjective optimization: A generalized homotopy approach_. Birkhäuser Basel, 2001.
  * [19] A. Martín and O. Schütze, “Pareto Tracer: a predictor–-corrector method for multi-objective optimization problems,” _Engineering Optimization_ , vol. 50, no. 3, pp. 516–536, 2018.
  * [20] O. Schütze, O. Cuate, A. Martín, S. Peitz, and M. Dellnitz, “Pareto Explorer: a global/local exploration tool for many-objective optimization problems,” _Engineering Optimization_ , vol. 52, no. 5, pp. 832–855, 2020.
  * [21] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv:1412.6980_ , 2014.
  * [22] A. Aşkan and S. Sayın, “SVM classification for imbalanced data sets using a multiobjective optimization framework,” _Annals of Operations Research_ , vol. 216, no. 1, pp. 191–203, 2014.
  * [23] H. Aytug and S. Sayın, “Exploring the trade-off between generalization and empirical errors in a one-norm SVM,” _European journal of operational research_ , vol. 218, no. 3, pp. 667–675, 2012.
  * [24] H. Li, X. Su, Z. Xu, and Q. Zhang, “MOEA/D with iterative thresholding algorithm for sparse optimization problems,” in _Parallel Problem Solving from Nature - PPSN XII_. Springer, 2012, pp. 93–101.
  * [25] M. Y. Park and T. Hastie, “L1-regularization path algorithm for generalized linear models,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 69, no. 4, pp. 659–677, 2007.
  * [26] S. Rosset, “Tracking curved regularized optimization solution paths,” in _Advances in Neural Information Processing Systems_. MIT Press, 2004.
  * [27] S. Rosset, “Following curved regularized optimization solution paths,” in _Advances in Neural Information Processing Systems_ , 2005, pp. 1153–1160.
  * [28] P. Zhao and B. Yu, “Boosted lasso,” University of California, Berkeley, Tech. Rep., 2004.
  * [29] F. H. Clarke, _Optimization and Nonsmooth Analysis_. Society for Industrial and Applied Mathematics, 1990.
  * [30] A. Bagirov, N. Karmitsa, and M. M. Mäkelä, _Introduction to Nonsmooth Optimization_. Springer International Publishing, 2014.
  * [31] M. M. Mäkelä, V. P. Eronen, and N. Karmitsa, “On nonsmooth multiobjective optimality conditions with generalized convexities,” in _Optimization in Science and Engineering: In Honor of the 60th Birthday of Panos M. Pardalos_ , 2014.
  * [32] E. L. Allgower and K. Georg, “Numerical continuation methods : An introduction,” 1990.
  * [33] J. Nocedal and S. Wright, _Numerical Optimization_. Springer New York, 2006.
  * [34] B. Gebken, S. Peitz, and M. Dellnitz, “On the hierarchical structure of Pareto critical sets,” _Journal of Global Optimization_ , vol. 73, no. 4, pp. 891––913, 2019.
  * [35] O. Schütze, “A New Data Structure for the Nondominance Problem in Multi-objective Optimization,” in _Lecture Notes in Computer Science_. Springer Berlin Heidelberg, 2003, pp. 509–518.
  * [36] R. A. Fisher, “The use of multiple measurements in taxonomic problems,” _Annals of Eugenics_ , vol. 7, no. 2, pp. 179–188, 1936.
  * [37] X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , 2010, pp. 249–256.
  * [38] Y. LeCun, C. Cortes, and C. Burges, “Mnist handwritten digit database,” _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , vol. 2, 2010.
  * [39] [Online]. Available: https://www.tensorflow.org/quantum/tutorials/mnist
  * [40] B. Gebken and S. Peitz, “An Efficient Descent Method for Locally Lipschitz Multiobjective Optimization Problems,” _Journal of Optimization Theory and Applications_ , vol. 188, pp. 696–723, 2021.
  * [41] J.-C. Loiseau and S. L. Brunton, “Constrained Sparse Galerkin Regression,” _Journal of Fluid Mechanics_ , vol. 838, pp. 42–67, 2018.

|  Katharina Bieker received the B.Sc. and M.Sc. degree in technomathematics
from Paderborn University, Germany, in 2015 and 2017, respectively. She is
currently a research assistant at the Chair of Applied Mathematics at
Paderborn University and is studying towards a PhD degree with a main focus on
data-driven modelling and multiobjective optimization.
---|---
|  Bennet Gebken received the B.Sc. and M.Sc. degree in mathematics from
Paderborn University, Germany, in 2014 and 2017, respectively. He is currently
a research assistant at the Chair of Applied Mathematics at Paderborn
University and is studying towards a PhD degree with a main focus on the
structure and computation of Pareto critical sets in multiobjective
optimization.
---|---
|  Sebastian Peitz received his B.Sc. and M.Sc. degrees in mechanical
engineering from RWTH Aachen University, Germany, in 2011 and 2013,
respectively, and his PhD in Mathematics from Paderborn University in 2017. He
is currently assistant professor for “Data Science for Engineering” in the
Department of Computer Science at Paderborn University. His research interests
are multiobjective optimization, optimal control and data-driven modeling of
complex systems.
---|---