[1]Michael Muehlebach

[1]Learning and Dynamical Systems, Max Planck Institute for Intelligent
Systems, Max-Planck-Ring 4, Tuebingen, 72076, Baden-Wuerttemberg, Germany

2]Department of Electrical Engineering and Computer Science, University of
California, Berkeley, 387 Soda Hall, Berkeley, 94720, California, USA

# Accelerated First-Order Optimization under Nonlinear Constraints

<EMAIL_ADDRESS>Michael I. Jordan<EMAIL_ADDRESS>*  [

###### Abstract

We exploit analogies between first-order algorithms for constrained
optimization and non-smooth dynamical systems to design a new class of
accelerated first-order algorithms for constrained optimization. Unlike Frank-
Wolfe or projected gradients, these algorithms avoid optimization over the
entire feasible set at each iteration. We prove convergence to stationary
points even in a nonconvex setting and we derive accelerated rates for the
convex setting both in continuous time, as well as in discrete time. An
important property of these algorithms is that constraints are expressed in
terms of velocities instead of positions, which naturally leads to sparse,
local and convex approximations of the feasible set (even if the feasible set
is nonconvex). Thus, the complexity tends to grow mildly in the number of
decision variables and in the number of constraints, which makes the
algorithms suitable for machine learning applications. We apply our algorithms
to a compressed sensing and a sparse regression problem, showing that we can
treat nonconvex $\ell^{p}$ constraints ($p<1$) efficiently, while recovering
state-of-the-art performance for $p=1$.

###### keywords:

Constrained Optimization, Nonlinear Programming, Gradient-based Methods,
Machine Learning

##  1 Introduction

Our work is concerned with developing first-order algorithms for nonlinear
constrained optimization problems of the following form:

$\min_{x\in\mathbb{R}^{n}}f(x),\quad\text{s.t.}\quad g(x)\geq 0,$ (1)

where the function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ defines the
objective, the function $g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n_{\text{g}}}$
the constraints, and where $n$ and $n_{\text{g}}$ are positive integers. The
applications we have in mind include problems in statistics, machine learning,
and control theory, where $n$ and $n_{\text{g}}$ are typically on the order of
$10^{6}$.

We develop a new class of accelerated first-order algorithms for computing
stationary points of (1). These algorithms have four distinctive features: (i)
they rely on local linear approximations of the feasible set, thereby avoiding
projections or optimizations over the entire feasible set at every iteration,
(ii) they are conceptually easy to understand and easy to implement, (iii)
their iteration complexity for convex problems is dimension-
independent,111This contrasts with interior-point methods, for example, which
require $\mathcal{O}(\sqrt{n_{\text{g}}})$ Newton iterations to decrease the
value of the objective by a constant factor. and (iv) in many important cases
(even nonconvex ones) their per-iteration complexity scales roughly linearly
in the number of decision variables and the number of constraints.

An important aspect of our work is to lift position constraints to a velocity
level, which naturally results in a local linear approximation of the feasible
set. These approximations come in two variants, each providing a different
trade-off between the per-iteration complexity and the resulting convergence
guarantees. More precisely, in our algorithms the forward increments,
$(x_{k+1}-x_{k})/T$, where $T>0$ is the step size, will be constrained to the
set

$V_{\alpha}(x):=\\{v\in\mathbb{R}^{n}\leavevmode\nobreak\
|\leavevmode\nobreak\ \nabla g_{i}(x)^{\mathsf{\tiny{T}}}v+\alpha g_{i}(x)\geq
0,\leavevmode\nobreak\ \forall i\in I\\},$ (2)

where $I$ either has the form $I=[n_{\text{g}}]:=\\{1,2,\dots,n_{\text{g}}\\}$
or the form $I_{x}:=\\{\leavevmode\nobreak\ i\in\mathbb{Z}\leavevmode\nobreak\
|\leavevmode\nobreak\ g_{i}(x)\leq 0\\}$. The former version includes every
constraint at every iteration, while the latter version includes only
constraints that are violated at the current iterate $x$. Our convergence
results span both continuous-time and discrete-time models and are summarized
as follows (see the corresponding theorems for the precise statements):

  * •

Our discrete-time algorithms achieves an accelerated linear rate (iteration
complexity $\mathcal{O}(\sqrt{\kappa_{l}}\log(1/\epsilon))$) if $f$ is smooth
and strongly convex, $g$ is smooth and concave, and if a bound on
$\lambda^{*}$, an optimal multiplier of (1) is known. The condition number
$\kappa_{l}$ refers to the Lagrangian
$f(x)-{\lambda^{*}}^{\mathsf{\tiny{T}}}g(x)$ and $\epsilon$ is the tolerance.

  * •

The same discrete-time algorithms achieve an accelerated sublinear rate
(iteration complexity $\mathcal{O}(1/\sqrt{\epsilon})$) if no bound on
$\lambda^{*}$ is known.

  * •

We also develop discrete-time algorithms for the nonconvex setting that ensure
convergence to stationary points provided that $f$ and $g$ are smooth.

  * •

Our continuous-time algorithms achieve an accelerated linear rate (iteration
complexity $\mathcal{O}(\sqrt{\kappa}\log(1/\epsilon))$ if $f$ is smooth and
strongly convex, $g$ is smooth and concave, where $\kappa$ denotes the
condition number of $f$.

  * •

If $f$ fails to be strongly convex, our continuous-time algorithms still
achieve an accelerated sublinear rate (iteration complexity
$\mathcal{O}(1/\sqrt{\epsilon})$).

  * •

If $f$ fails to be convex and $g$ fails to be concave we show convergence to
the set of stationary points of our continuous-time algorithms.

In the special case where $g$ is linear, the condition number $\kappa_{l}$
reduces to $\kappa$ and our algorithm recovers the rate of accelerated
projected gradient descent, without requiring projections or optimizations
over the entire feasible set.

Our treatment builds on recent progress in using tools from continuous-time
dynamical systems to analyze discrete-time algorithms in gradient-based
optimization [1, 2, 3, 4, 5, 6, 7, 8, 9]. Much of this work aims at
understanding accelerated first-order optimization methods, such as Nesterov’s
algorithm, by exposing links between differential and symplectic geometry,
dynamical systems, and mechanics. While in the absence of constraints these
analogies result in _smooth_ dynamical systems, the current article presents
analogies between constrained optimization and _non-smooth_ dynamical systems.
Indeed, one of the closest point of contacts with existing literature is the
notion of Moreau time-stepping in non-smooth mechanics [10]. The important
feature of Moreau time-stepping, which also lies at the heart of our work, is
that smooth and non-smooth motion are treated on equal footing, which is
achieved by discretizing a certain kind of differential inclusion [see, e.g.,
11, 12].

Our approach can also be interpreted through the lens of the projected
gradient methodology and indeed it has certain similarities to inexact
projected gradient methods, as proposed by Wang and Liu [13] and Birgin et al.
[14]. While projected gradient approaches have been successfully applied in
various machine learning problems [see, e.g., 15, 16], the Frank-Wolfe
algorithm has also received considerable attention in recent years [17]. The
appeal of Frank-Wolfe is further increased by the fact that it provides a
unified framework for many first-order machine learning algorithms in
constrained settings, including support vector machines, online estimation of
mixtures of probability densities, and boosting [18]. Recent results extend
the Frank-Wolfe algorithm to the stochastic setting [19, 20], or improve on
its relatively slow convergence rate [21, 22].

In some cases constraints can be handled very efficiently with mirror descent,
[23, Ch. 3], where a non-Euclidean metric is introduced that adapts gradient
descent to the specific type of objective function or the specific type of
constraints at hand [24]. Although mirror descent is based on projections onto
the feasible set, the non-Euclidean metric can improve on problem-specific
constants. An important example is the optimization of linear functions over
the unit simplex, which has applications in online machine learning [25].

Compared to projected gradients, mirror descent, and the Frank-Wolfe
algorithm, our approach avoids optimizing over the entire feasible set at each
iteration and instead relies on sparse, local and convex approximations. This
article focuses on accelerated gradient descent, building on the recent
results of Muehlebach and Jordan [26], which treats gradient descent.
Including constraints in momentum-based algorithms is challenging: The
presence of constraints requires a need for sudden and large changes in
momentum (impacts) in order to avoid infeasible iterates. This requires us to
not only characterize the smooth motion (if constraints are absent or the
solution slides along the boundary of the feasible set), but also the non-
smooth motion (if the solution suddenly hits the boundary of the feasible
set).

Specific problems which have the potential to benefit from our approach
include planning problems in reinforcement learning and/or optimal control
[see, e.g., 27], optimizations over nonconvex matrix manifolds (such as the
set of orthogonal matrices [see, e.g., 28]), distance geometry problems that
arise in computational chemistry/NMR spectroscopy [29], optimal transport
problems [30], or supervised learning tasks that involve nonlinear constraints
(for example in an imitation learning framework, where nonlinear constraints
arise from stability requirements on the closed-loop system). We will also
demonstrate our approach on $\ell^{p}$-regularized inverse problems that arise
in compressed sensing and signal processing, where we are not only able to
obtain state-of-the-art results for $p=1$, but can also seamlessly handle the
regime $0<p<1$.

The article is structured in the following way: Sec. 2 summarizes earlier work
of Muehlebach and Jordan [26], which covers gradient descent and sets the
stage for discussing momentum-based algorithms in Sec. 3. A variety of
convergence results that capture both discrete-time and continuous-time models
are presented in Sec. 4; in particular, in the nonconvex regime we establish
convergence to stationary points and we derive accelerated rates in the convex
regime. Sec. 5 presents numerical experiments, which include nonconvex sparse
regression and compressed sensing problems. The paper concludes with a short
discussion in Sec. 6.

##  2 Constrained Gradient Flow

The fundamental idea of this work is to express constraints in terms of the
forward increment or velocity of our algorithms instead of constraining the
iterates or positions directly. As we will see shortly, this naturally leads
to local, sparse and convex approximations of the feasible set. Our treatment
builds upon Muehlebach and Jordan [26], which focused on gradient descent and
gradient flow, whereas this work focuses on accelerated gradient algorithms.
In order to simplify the presentation, we begin with a brief review of this
earlier work. The section ends with deriving a non-asymptotic linear
convergence rate for the resulting gradient-descent algorithm, something which
was missing in [26]. This will also set the stage for deriving accelerated
rates in Sec. 4.

We model an optimization algorithm as a continuous-time or discrete-time
dynamical system, whose equilibria correspond to the stationary points of (1).
In continuous time, the configuration of the system will be denoted by a
function $x:[0,\infty)\rightarrow\mathbb{R}^{n}$, which is assumed to be
absolutely continuous. A fundamental observation, lying at the heart of the
current research, is that the constraint $x(t)\in C$, for all $t\geq 0$, is
equivalent to the constraint $\dot{x}(t)^{+}\in T_{C}(x(t))$, for all $t\geq
0$, $x(0)\in C$, where $T_{C}(x(t))$ denotes the tangent cone (in the sense of
Clarke) of the set $C$ at $x(t)\in\mathbb{R}^{n}$, and $\dot{x}(t)^{+}$
denotes the forward velocity: $\dot{x}(t)^{+}:=\lim_{\mathrm{d}t\downarrow
0}(x(t+\mathrm{d}t)-x(t))/\mathrm{d}t$. The tangent cone $T_{C}(x)$ is defined
as the set of all vectors $v$ such that $(x_{k}-x)/t_{k}\rightarrow v$ for two
sequences $x_{k}\in C$ and $t_{k}\geq 0$ with $x_{k}\rightarrow x$,
$t_{k}\rightarrow 0$. Provided that a constraint qualification holds (for
example Mangasarian-Fromovitz or Abadie constraint qualification), the tangent
cone can be expressed as

$T_{C}(x)=\\{v\in\mathbb{R}^{n}\leavevmode\nobreak\ |\leavevmode\nobreak\
\nabla g_{i}(x)^{\mathsf{\tiny{T}}}v\geq 0,\leavevmode\nobreak\
\leavevmode\nobreak\ \forall i\in I_{x}\\},$

where $I_{x}=\\{i\in\mathbb{Z}\leavevmode\nobreak\ |\leavevmode\nobreak\
g_{i}(x)\leq 0\\}$ denotes the set of active inequality constraints at $x$.

We therefore conclude that the constraint $x(t)\in C$, which constrains the
position $x$, is equivalent to a constraint on the forward velocity
$\dot{x}^{+}$. We note that the velocity $\dot{x}$ is allowed to be
discontinuous and need not exist for every $t\geq 0$.222We assume that
$\dot{x}$ is of locally bounded variation, which means that on any compact
interval $\dot{x}$ has countably many discontinuity points, where left and
right limits exist. For example, if the trajectory $x$ reaches the boundary of
the feasible set, an instantaneous jump of the velocity might be required to
ensure that $x$ remains in $C$.

In discrete time, however, this equivalence between position and velocity
constraints no longer holds, since $T_{C}(x)$ is only a first-order
approximation of the feasible set. Thus, implementing $(x_{k+1}-x_{k})/T\in
T_{C}(x_{k})$ may lead to infeasible iterates. Muehlebach and Jordan [26]
therefore suggest to introduce the velocity constraint $V_{\alpha}(x)$, see
(2) with $I=I_{x}$, which includes the restitution coefficient $\alpha>0$. The
following remarks motivate (2):

  * i)

For $x\in C$, the set $V_{\alpha}(x)$ reduces to the tangent cone $T_{C}(x)$
(assuming constraint qualification).

  * ii)

For a fixed $x\in\mathbb{R}^{n}$, $V_{\alpha}(x)$ is a convex polyhedral set
involving only the active constraints $I_{x}$. The set $V_{\alpha}(x)$
therefore amounts to a sparse and linear approximation of the feasible set
$C$, even if $C$ is nonconvex.

  * iii)

In continuous time, the constraint $\dot{x}(t)^{+}\in V_{\alpha}(x(t))$ for
all $t\geq 0$ implies

$\displaystyle g_{i}(x(t))$ $\displaystyle\geq\min\\{g_{i}(x(0))e^{-\alpha
t},0\\},$ (3)

for all $t\geq 0$ and all $i\in\\{1,\dots,n_{\text{g}}\\}$, which can be
verified with Grönwall’s inequality. This means that potential constraint
violations decrease at rate $\alpha$.

The continuous-time gradient flow dynamics that were studied in Muehlebach and
Jordan [26] arise from the following equation:

$\dot{x}(t)^{+}+\nabla f(x(t))=R(t),\quad-R(t)\in
N_{V_{\alpha}(x(t))}(\dot{x}(t)^{+}),$ (4)

for all $t\geq 0$, where $N_{V_{\alpha}(x(t))}(\dot{x}(t)^{+})$ denotes the
normal cone of the set $V_{\alpha}(x(t))$ at $\dot{x}(t)^{+}$. Thus, the
variable $R(t)$ can be regarded as a constraint force that imposes the
constraint $\dot{x}(t)^{+}\in V_{\alpha}(x(t))$. Moreover, for each fixed
$t\geq 0$, we can eliminate $R(t)$ in (4) and interpret the resulting
expression as a stationarity condition with respect to $\dot{x}(t)^{+}$, which
means that (4) is equivalent to

$\dot{x}(t)^{+}:=\mathop{\mathrm{argmin}}_{v\in
V_{\alpha}(x(t))}\frac{1}{2}|v+\nabla f(x(t))|^{2}.$ (5)

As long as $V_{\alpha}(x(t))$ is nonempty, this guarantees the uniqueness of
$\dot{x}(t)^{+}$ for every $t\geq 0$. It also provides the following natural
interpretation: the forward velocity $\dot{x}(t)^{+}$ is chosen to match the
unconstrained gradient flow equation $\dot{x}(t)^{+}+\nabla f(x(t))=0$ as
nearly as possible, subject to the velocity constraint $\dot{x}(t)^{+}\in
V_{\alpha}(x(t))$.

######  Remark 1.

Nonemptiness of $V_{\alpha}(x)$: If $C$ is convex, $V_{\alpha}(x)$ is
guaranteed to be nonempty for all $x\in\mathbb{R}^{n}$. If $C$ is nonconvex,
nonemptiness of $V_{\alpha}(x)$ for all $x$ in a neighborhood of $C$ is
guaranteed if the Mangasarian-Fromovitz constraint qualification holds for all
$x\in C$. We note that the Mangasarian-Fromovitz constraint qualification is
generic in the following sense: Provided that $g$ is semi-algebraic (these
cases include all the usual functions used in optimization) there exists
$\epsilon_{0}\in\mathbb{R}^{m},\epsilon_{0}>0$, such that the set
$C_{\epsilon}:=\\{x\in\mathbb{R}^{n}\leavevmode\nobreak\ |\leavevmode\nobreak\
g(x)\geq-\epsilon\\}$ satisfies the Mangasarian-Fromovitz constraint
qualification for all $x\in C_{\epsilon}$ and for all
$\epsilon\in(0,\epsilon_{0})$; see Bolte et al. [31].

In discrete time, it suffices to replace $\dot{x}(t)^{+}$ by
$(x_{k+1}-x_{k})/T$ and $x(t)$ by $x_{k}$ in order to obtain the corresponding
constrained gradient-descent dynamics. The condition (5) can be interpreted as
a modified projected gradient scheme, where projections over the entire
feasible set $C$ are replaced with optimizations over the sparse and convex
approximation $V_{\alpha}(x_{k})$. The remark about the nonemptiness of
$V_{\alpha}(x_{k})$ applies in the same way.

The results from Muehlebach and Jordan [26] establish convergence of (4)
(and/or (5)) both in continuous and discrete time. In continuous time, it was
shown that even when $f$ and $C$ are nonconvex, the trajectories of (4)
(and/or (5)) converge to the set of stationary points. Moreover, if $f$ is
strongly convex with strong convexity constant $\mu$ and $\alpha$ is set to
$2\mu$, the trajectories converge from any initial condition to the minimizer
of (1) at a linear rate:

$|f(x(t))-f(x^{*})|\leq(|f(x(0))-f^{*}|+c_{1})e^{-2t/\kappa},$

where $\kappa$ is the condition number of $f$ and $c_{1}\geq 0$ is an explicit
constant that captures whether the initial condition is feasible or not. The
notation $|\cdot|$ is used throughout the article to denote the $\ell^{2}$
norm. Similarly, if $f$ is strongly convex and $C$ is convex, in discrete time
(where $\dot{x}(t)^{+}$, $x(t)$, and $R(t)$ are replaced with
$(x_{k+1}-x_{k})/T$, $x_{k}$, and $R_{k}$, respectively) the trajectories were
shown to converge to the minimizer of (1) for $T\leq 2/(L_{\text{l}}+\mu)$,
$\alpha<\mu$, where $L_{\text{l}}$ refers to the smoothness constant of the
corresponding Lagrangian. The asymptotic convergence was found to be linear at
a rate of $1/\kappa_{l}$, where $\kappa_{l}$ denotes to the Lipschitz constant
of the Lagrangian. However, a non-asymptotic linear convergence result was
missing in discrete time; the asymptotic result required a relatively complex
analysis of the algorithm’s dynamics. We now provide such an non-asymptotic
result.

######  Theorem 1.

Let the function $f$ be $\mu$-strongly convex and $g$ be smooth and concave.
Let $l(x):=f(x)-{\lambda^{*}}^{\mathsf{\tiny{T}}}g(x)$ denote the Lagrangian,
where $\lambda^{*}$ denotes an optimal multiplier of (1), and let $T\leq
1/L_{l}$, $\alpha\leq\mu$. Then, the iterates of Alg. 1 satisfy:

$\frac{\mu}{2}|x_{k}-x^{*}|^{2}\leavevmode\nobreak\ \leavevmode\nobreak\
\leq\leavevmode\nobreak\ \leavevmode\nobreak\
l(x_{k})-l(x^{*})\leavevmode\nobreak\ \leavevmode\nobreak\
\leq\leavevmode\nobreak\ \leavevmode\nobreak\
(1-T\alpha)^{k}\leavevmode\nobreak\ \leavevmode\nobreak\ (l(x_{0})-l(x^{*})),$

for all $k\geq 0$, where $x^{*}$ denotes the minimizer of (1). For $T=1/L_{l}$
and $\alpha=\mu$, the convergence is linear at rate $\kappa_{l}=L_{l}/\mu$,
where $\mu$ and $L_{l}$ denote the strong convexity and smoothness constants
of $l$.

Algorithm 1 Constrained gradient descent with velocity projections

$x_{0}\in\mathbb{R}^{n}$, $L_{l},\mu$ (smoothness and strong convexity
constant of $l$)

MAXITER (maximum number of iterations)

$T\leq 1/L_{l}$, $\alpha\leq\mu$

for $k=0,1,\dots$ do

$x_{k+1}\leftarrow\mathop{\mathrm{argmin}}_{v\in V_{\alpha}(x_{k})}|v+\nabla
f(x_{k})|^{2}$ (with $I=[n_{\text{g}}]$)

end for

The proof of Thm. 1 is included in App. A, where we also add a corollary that
shows that a $1/k$ rate is obtained if the smoothness constant $L_{l}$ is
unknown (knowledge of $L_{l}$ requires a bound on $\lambda^{*}$).

This implies that (4) (and/or (5)) implement gradient-flow and gradient-
descent dynamics that can handle constraints and converge linearly with the
typical $1/\kappa$ and $1/\kappa_{l}$-rate if the objective function $f$ is
smooth and strongly convex and a bound on the optimal multiplier $\lambda^{*}$
is known. The set $V_{\alpha}(x)$ can be seen as a velocity constraint and
provides a natural generalization of the tangent cone. It also reduces the
computational cost for each iteration, since projections on the entire
feasible set are avoided. In the next section, we generalize these ideas to
algorithms that have momentum. This will naturally lead to accelerated
algorithms that converge linearly at a rate of $1/\sqrt{\kappa_{l}}$ (if $f$
is smooth and strongly convex) or at the sublinear rate $1/t^{2}$ (if $f$ is
smooth and convex), which is a significant speedup. We will also derive
convergence if $f$ and $C$ are nonconvex.

##  3 Accelerated Gradient Flow

We begin our presentation with a derivation in continuous time. The
corresponding discrete-time algorithms will be stated subsequently. A natural
starting point is the work of authors such as Polyak [32], Su et al. [1], and
Muehlebach and Jordan [7], who argued that in the _unconstrained_ case,
accelerated optimization algorithms can be viewed as dynamical systems
described by second-order differential equations. A canonical example is the
following:

$\dot{u}(t)+2\delta u(t)+\nabla f(x(t)+\beta u(t))=0,$ (6)

where we use the variable $u(t)=\dot{x}(t)$ to denote the velocity (or
momentum), and where $\delta\geq 0$ and $\beta\geq 0$ are damping
parameters.333The variables $\delta$, $\beta$ may also depend on time. For
ease of presentation we focus on the case where $\delta$ and $\beta$ are
fixed, but also state corresponding results for time-varying parameters.

In the presence of constraints, $u(t)$ is allowed to be discontinuous, which
is in sharp contrast to (6). For example, if the trajectory $x(t)$ approaches
the boundary of the feasible set, an instantaneous jump in $u(t)$ might be
required to ensure that $x(t)$ remains feasible. Thus, compared to (4) (or
equivalently (5)), where the state $x(t)$ is absolutely continuous, we are now
in a position where we allow for the state $(x(t),u(t))$ (which includes the
velocity $u$) to be discontinuous. This means that in addition to a
differential equation of the type (4), which characterizes the smooth motion,
we also prescribe how the discontinuities in $u$ can arise. If we regard
$(x(t),u(t))$ as the position and velocity of a mechanical system,
discontinuities in $u$ have a mechanical meaning as impacts, which are
described by a corresponding impact law. The mathematical formalism, which
enables discontinuities in $u$, is summarized next.

We regard the state $z:=(x,u)$ to be the result of an integration process:

$z(t)=z(t_{0})+\int_{t_{0}}^{t}\mathrm{d}z,\quad\forall t\geq t_{0}.$

However, instead of the usual Lebesgue density
$\mathrm{d}z=\dot{z}(t)\mathrm{d}t$, $\mathrm{d}z$ now represents a
differential measure [33], and admits both a density with respect to the
Lebesgue measure (denoted by $\mathrm{d}t$), as well as a density with respect
to an atomic measure (denoted by $\mathrm{d}\eta$). As is common in non-smooth
mechanics, we assume that $z(t)$ is of locally bounded variation and does not
contain any singular terms. This means that $z(t)$ can be decomposed as an
absolutely continuous function and a piecewise constant step function [33]. At
every time $t$, $z(t)$ has well-defined left and right limits, $z(t)^{-}$ and
$z(t)^{+}$, even though $z(t)$ might not exist or might not be of interest. We
can express the differential measure $\mathrm{d}z$ as
$\mathrm{d}z=\dot{z}(t)\mathrm{d}t+(z(t)^{+}-z(t)^{-})\mathrm{d}\eta,$ and the
integration over an interval $[t_{0},t]$, which contains the time instants
$t_{\text{d}i}$, $i=1,2,\dots$, where $z(t)$ is discontinuous, yields

$z(t)^{+}=z(t_{0})^{-}+\int_{t_{0}}^{t}\dot{z}(t)\mathrm{d}t+\sum_{i\geq
1}z(t_{\text{d}i})^{+}-z(t_{\text{d}i})^{-}.$

As a consequence of allowing the state to be discontinuous, we need to
delineate both the density $\dot{z}(t)$ with respect to the Lebesgue measure
$\mathrm{d}t$ (which describes the smooth part of the motion) as well as the
density $z(t)^{+}-z(t)^{-}$ (which describes the non-smooth motion) for fully
determining the state trajectory $z(t)$. By analogy to non-smooth mechanics
[see, e.g., 12], this can be achieved with the following measure-differential
inclusion:

$\mathrm{d}u+2\delta u\mathrm{d}t+\nabla f(x+\beta u)\mathrm{d}t=\sum_{i\in
I_{x}}\nabla
g_{i}(x)\mathrm{d}\lambda_{i},\quad\gamma_{i}^{+}+\epsilon\gamma_{i}^{-}\in
N_{\mathbb{R}_{\leq 0}}(-\mathrm{d}\lambda_{i}),\quad i\in I_{x},$ (7)

where $\epsilon\in[0,1)$ is a constant, $\gamma_{i}$ is the velocity
associated with the $i$th constraint and is defined as

$\gamma_{i}(x,u):=\nabla g_{i}(x)^{\mathsf{\tiny{T}}}u+\alpha g_{i}(x),$

and where we have omitted the dependence on $t$ (as we will do frequently in
the subsequent presentation). We note that the set $I_{x}$ (or $I_{x(t)}$ in
full notation) is time-dependent. The normal cone inclusion in (7) is
illustrated with Fig. 1 and will be further discussed below. The constant
$\epsilon$ has the interpretation of a restitution coefficient, whereby
$\epsilon=0$ leads to inelastic collisions, and $\epsilon=1$ yields elastic
collisions. Measure-differential inclusions are common in non-smooth
mechanics, where various existence results for inclusions of the type (7) have
been established; see, for example, Piazza et al. [34], Leine and van de Wouw
[35] and references therein.

We note that if $x$ is in the interior of the feasible set, $I_{x}$ is empty,
and therefore (7) reduces to (6). This means that (7) generalizes (6) from the
unconstrained case to the constrained case by including the constraint
$\gamma_{i}^{+}+\epsilon\gamma_{i}^{-}\in N_{\mathbb{R}_{\leq
0}}(-\mathrm{d}\lambda_{i})$, which, as we will discuss below, describes the
discontinuities of $u$ via Newton’s impact law and imposes the velocity
constraint $u(t)\in V_{\alpha}(x(t))$, whenever $u(t)$ exists.

$\gamma_{i}^{+}+\epsilon_{i}\gamma_{i}^{-}$$-\mathrm{d}\lambda_{i}$

$-0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.5$$1$$x$$(x)_{\Delta}^{p},\leavevmode\nobreak\
\leavevmode\nobreak\ x^{p}$$(x)_{\Delta}^{p}$$x^{p}$

Figure 1: The left panel shows the normal cone inclusion
$\gamma_{i}^{+}+\epsilon\gamma_{i}^{-}\in N_{\mathbb{R}_{\leq
0}}(-\mathrm{d}\lambda_{i})$, which is equivalent to the complementarity
condition $\mathrm{d}\lambda_{i}\geq 0$,
$\gamma_{i}^{+}+\epsilon\gamma_{i}^{-}\geq 0$,
$\mathrm{d}\lambda_{i}(\gamma_{i}^{+}+\epsilon\gamma_{i}^{-})=0$. The right
panel shows the approximation $(x)^{p}_{\Delta}$ of $x^{p}$ for $\Delta=0.01$
and $p=0.6$. There is an excellent agreement between the approximation and
$x^{p}$ even though $\Delta$ is comparably large. In the numerical
experiments, see Sec. 5, $\Delta$ is set to $10^{-6}$.

It is important to note that (7) is understood in the sense of integration:
For any compact time interval $[t_{0},t_{1}]$, (7) defines the difference
$u(t_{1})^{+}-u(t_{0})^{-}$, which is obtained by integrating $\mathrm{d}u$
from $t_{0}$ to $t_{1}$; similarly, the difference $x(t_{1})^{+}-x(t_{0})^{-}$
is obtained by integrating $u(t)\mathrm{d}t$. This means that (7) has a very
natural discretization, which will be discussed in the next paragraph. The
last two paragraphs describe the interpretation of (7) in terms of its smooth
and non-smooth components. Formal convergence results in continuous and
discrete time will be derived in Sec. 4.

#### Discretization of (7):

The measure-differential inclusion (7) lends itself to the following
discretization: $\mathrm{d}u=u_{k+1}-u_{k}$, $\mathrm{d}t=T_{k}$,
$\mathrm{d}\lambda_{i}=\Lambda_{ki}$,
$\gamma_{i}^{+}=\gamma_{i}(x_{k},u_{k+1})$,
$\gamma_{i}^{-}=\min\\{0,\gamma_{i}(x_{k},u_{k})\\}$, where $T_{k}>0$ is the
step size.444The $\min$ in $\gamma_{i}^{-}$ ensures $u_{k+1}\in
V_{\alpha}(x_{k})$, which is not automatically satisfied in discrete time.
This yields

$\displaystyle u_{k+1}-u_{k}+\delta u_{k}T_{k}+\nabla f(x_{k}+\beta
u_{k})T_{k}$ $\displaystyle=\sum_{i\in I_{x_{k}}}\nabla
g_{i}(x_{k})\Lambda_{ki},$
$\displaystyle\gamma_{i}(x_{k},u_{k+1})+\epsilon\min\\{0,\gamma_{i}(x_{k},u_{k})\\}$
$\displaystyle\in N_{\mathbb{R}_{\leq 0}}(-\Lambda_{ki}),\quad i\in
I_{x_{k}}.$

We use the newly computed momentum for updating the position $x_{k}$:
$x_{k+1}=x_{k}+T_{k}u_{k+1}$, which is motivated by analogy to unconstrained
optimization. (This discretization scheme is found to be superior compared to
the standard Euler method [9].) The resulting update for $u_{k+1}$ can be
interpreted as a stationarity condition for $u_{k+1}$, and as a result, the
proposed algorithm can be summarized as follows:

$\displaystyle u_{k+1}$
$\displaystyle\\!=\\!\mathop{\mathrm{argmin}}_{v\in\mathbb{R}^{n}}\\!\frac{1}{2}|v\\!-\\!u_{k}\\!+\\!2\delta
u_{k}T_{k}\\!+\\!\nabla f(x_{k}+\beta u_{k})T_{k}|^{2},$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\text{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
\gamma_{i}(x_{k},v)\geq-\epsilon\min\\{\gamma_{i}(x_{k},u_{k}),0\\},\leavevmode\nobreak\
\leavevmode\nobreak\ i\in I_{x_{k}}$ $\displaystyle x_{k+1}$
$\displaystyle\\!=\\!x_{k}+T_{k}u_{k+1}.$ (8)

Remark 1 applies here in the same way: If $C$ is convex and $\epsilon=0$, the
feasible set in (8) is guaranteed to be nonempty, which means that $u_{k+1}$
is well defined (existence and uniqueness). If $C$ is nonconvex or
$\epsilon>0$, nonemptiness of the feasible set is guaranteed if constraint
qualifications are satisfied (for example Mangasarian-Fromovitz). These
constraint qualifications are generic, as discussed in Remark 1, which ensures
that $u_{k+1}$ is well defined as long as $x_{k}$ stays in a neighborhood of
the feasible set. As will be shown with our convergence analysis (see Sec. 4),
we can indeed ensure that $x_{k}$ remains in a neighborhood of $C$; the size
of the neighborhood can be controlled by the value of $T_{k}$. The pseudo-code
of the full algorithm is listed in App. B.

The following remarks are important:

  * i)

The update (8) has the interpretation of choosing $u_{k+1}$ to be as close as
possible to the update in the unconstrained case subject to the velocity
constraint
$\gamma_{i}(x_{k},u_{k+1})\geq-\epsilon\min\\{\gamma_{i}(x_{k},u_{k}),0\\}$.
As a result, in case $I_{x_{k}}$ is empty, (8) reduces to a standard momentum-
based method; if $\beta=0$ we obtain the heavy-ball algorithm, if $\beta\neq
0$ we obtain Nesterov’s method.

  * ii)

The update (8) includes only the constraints $I_{x_{k}}$ which are active at
iteration $k$. The constraint on $v$ in (8) is guaranteed to be convex, even
if the underlying feasible set is nonconvex. The constraints in (8) yield
therefore a sparse, local and convex approximation of the feasible set.
Instead of performing optimizations on the position level as is common with
projected gradients or the Frank-Wolfe method, (8) suggests to constrain the
velocities $u_{k}$, $k=1,2,\dots$.

We now proceed to give an interpretation and explanation of the continuous-
time dynamics (7).

#### Smooth motion:

If $u(t)$ happens to be absolutely continuous in the interval $(t_{0},t_{1})$,
its differential measure reduces to $\dot{u}(t)\mathrm{d}t$. Similarly, the
multipliers $\mathrm{d}\lambda_{i}$ have only a density with respect to the
Lebesgue measure $\mathrm{d}t$, which we denote by $\lambda_{i}(t)$, i.e.,
$\mathrm{d}\lambda_{i}=\lambda_{i}(t)\mathrm{d}t$. As a result, (7) reduces to

$\dot{u}+2\delta u+\nabla f(x+\beta u)=\sum_{i\in I_{x}}\nabla
g_{i}(x)\lambda_{i},$ (9)

for all $t\in(t_{0},t_{1})$ (a.e.). Furthermore, absolute continuity of $u(t)$
implies absolute continuity of $\gamma_{i}$, i.e.,
$\gamma_{i}^{+}=\gamma_{i}^{-}=\gamma_{i}$. In the limit
$\mathrm{d}t\downarrow 0$, the inclusion in (7) therefore reduces to

$(1+\epsilon)\gamma_{i}\in N_{\mathbb{R}_{\leq
0}}(-\lambda_{i})\quad\Leftrightarrow\quad\gamma_{i}\in N_{\mathbb{R}_{\leq
0}}(-\lambda_{i}),$

for all $i\in I_{x(t)}$ and for all $t\in(t_{0},t_{1})$ (a.e.).
($N_{\mathbb{R}_{\leq 0}}$ is a cone; we can therefore divide by
$1+\epsilon>0$.) The normal cone inclusion prescribing the relationship
between $\gamma_{i}$ and $\lambda_{i}$ is similar to Fig. 1.

From a physics perspective the normal cone inclusion $\gamma_{i}\in
N_{\mathbb{R}\leq 0}(-\lambda_{i})$ represents a force law, which by conic
duality can also be expressed as (see again Fig. 1)

$-\lambda_{i}\in N_{\mathbb{R}\geq 0}(\gamma_{i}).$

The sum $\nabla g_{i}(x(t))\lambda_{i}$ over $i\in I_{x}$ on the right-hand
side of (9) therefore has a physical interpretation as a constraint force:

$-R\\!=\\!-\\!\sum_{i\in I_{x}}\\!\nabla g_{i}(x)\lambda_{i},\quad\\!-R\in
N_{V_{\alpha}(x)}(u),$

which imposes the velocity constraint $u(t)\in V_{\alpha}(x(t))$ for all
$t\in(t_{0},t_{1})$ (a.e.). By virtue of Grönwall’s inequality this ensures
that constraint violations decrease linearly at rate $\alpha$, as highlighted
in (3).

We therefore conclude that in case of smooth motion, the measure-differential
inclusion (7) generalizes the differential equation (6) from the unconstrained
case to the constrained case, where the additional constraint force $R(t)$
imposes the velocity constraint $u(t)\in V_{\alpha}(x(t))$ (for almost all
$t$). The introduction of the force $R(t)$ is analogous to (4).

Since the motion is smooth for almost every $t$, the normal cone inclusion in
(7) guarantees the satisfaction of the velocity constraint $u(t)\in
V_{\alpha}(x(t))$ (or equivalently, $\gamma_{i}(x(t),u(t))\geq 0$ for all
$i\in I_{x(t)}$) for all $t\geq 0$ (a.e.). However, when a new constraint
arises at time $t_{0}$, there might be a situation where
$\gamma_{i}^{-}(x(t_{0}),u(t_{0}))<0$. In such a case an impact will be
required to ensure that $\gamma_{i}^{+}(x(t_{0}),u(t_{0}))\geq 0$. This is the
subject of the next paragraph.

#### Non-smooth motion:

In order to derive the non-smooth motion we integrate (7) over a time instant
$\\{t\\}$, where $u(t)$ is discontinuous; that is, $u(t)^{-}\neq u(t)^{+}$.
Due to the fact that the singleton $\\{t\\}$ has zero Lebesgue measure, we are
left with the atomic parts, leading to
$\mathrm{d}u=(u(t)^{+}-u(t)^{-})\mathrm{d}\eta$,
$\mathrm{d}\lambda_{i}:=\Lambda_{i}\mathrm{d}\eta$,

$u(t)^{+}-u(t)^{-}=\sum_{i\in I_{x(t)}}\nabla
g_{i}(x(t))\Lambda_{i},\quad\gamma_{i}^{+}+\epsilon\gamma_{i}^{-}\in
N_{\mathbb{R}\leq 0}(-\Lambda_{i}),\quad i\in I_{x(t)}.$ (10)

The normal cone inclusion should be interpreted as a generalization of
Newton’s impact law. For $\Lambda_{i}>0$, it implies
$\gamma_{i}^{+}+\epsilon\gamma_{i}^{-}=0$, meaning that the velocity
associated to constraint $i$ after impact, $\gamma_{i}^{+}$, is
$-\epsilon\gamma_{i}^{-}$, where $\gamma_{i}^{-}$ is the velocity associated
to constraint $i$ before impact. From the discussion of the smooth motion it
follows $\gamma_{i}(x(t_{0}),u(t_{0}))^{-}$ at time $t_{0}$ can only be
negative if the constraint $i$ becomes active at time $t_{0}$; that is,
$i\not\in I_{x(t)}$ for $t<t_{0}$ and $i\in I_{x(t)}$ for $t=t_{0}$. This
necessitates a discontinuity in $u$ at time $t_{0}$, which according to the
above normal cone inclusion comes in two variants: i) $\Lambda_{i}>0$, which
implies $\gamma_{i}^{+}=-\epsilon\gamma_{i}^{-}$ and ii) $\Lambda_{i}=0$,
which implies $\gamma_{i}^{+}\geq-\epsilon\gamma_{i}^{-}$. In variant i), the
impulsive force $\Lambda_{i}$ contributes the component $\Lambda_{i}\nabla
g_{i}(x(t))$ (normal to constraint $i$) to the velocity jump
$u(t_{0})^{+}-u(t_{0})^{-}$, whereas in variant ii), there is no such
contribution. Both variants ensure
$\gamma_{i}(x(t_{0}),u(t_{0}))^{+}\geq-\epsilon\gamma_{i}(x(t_{0}),u(t_{0}))^{-}\geq
0$.

The characterization of the non-smooth motion according to (10) can be
interpreted as a stationarity condition for $u(t)^{+}$, which yields

$u(t)^{+}=\mathop{\mathrm{argmin}}_{v\in\mathbb{R}^{n}}\frac{1}{2}|v-u(t)^{-}|^{2}\quad\text{s.t.}\leavevmode\nobreak\
\leavevmode\nobreak\
\gamma_{i}(x(t),v)\geq-\epsilon\gamma_{i}(x(t),u(t)^{-}),\leavevmode\nobreak\
\leavevmode\nobreak\ \forall i\in I_{x(t)}.$ (11)

The minimization in (11) has the following meaning: for each $u(t)^{-}$ there
is a unique $u(t)^{+}$, which is chosen to be as close as possible to
$u(t)^{-}$ subject to Newton’s impact law
$\gamma_{i}^{+}\geq-\epsilon\gamma_{i}^{-}$ for all $i\in I_{x(t)}$.

#### Equilibria of (7):

The equilibria of (7) are obtained from $x(t)\equiv x_{0}$,
$\mathrm{d}\lambda_{i}\equiv\lambda_{0i}\mathrm{d}t$, $u(t)\equiv 0$,
$\mathrm{d}u\equiv 0$, where $x_{0}\in\mathbb{R}^{n}$ and the multipliers
$\lambda_{0i}\geq 0$, $i\in I_{x_{0}}$ are constant. As a result, (7) reduces
to

$-\nabla f(x_{0})+\sum_{i\in I_{x_{0}}}\nabla
g_{i}(x_{0})\lambda_{i0}=0,\leavevmode\nobreak\ \leavevmode\nobreak\
(1+\epsilon)\alpha g_{i}(x_{0})\in N_{\mathbb{R}_{\leq
0}}(-\lambda_{i0}),\leavevmode\nobreak\ \leavevmode\nobreak\ i\in I_{x_{0}}.$

The normal cone inclusion can be simplified by dividing by
$\alpha(1+\epsilon)>0$ (the normal cone is a cone), which implies that
$g_{i}(x_{0})$ and $\lambda_{i0}$ satisfy the complementarity conditions

$g_{i}(x_{0})\geq 0,\quad\lambda_{i0}\geq
0,\quad\lambda_{i0}g_{i}(x_{0})=0,\quad\forall i\in I_{x_{0}}.$

Hence, the equilibria of (7) satisfy the Karush-Kuhn-Tucker conditions of (1),
which means that the stationary points of (1) are indeed equilibria.

##  4 Convergence Analysis

The following section discusses the convergence of trajectories of (7) and
(8), and characterizes the rate of convergence. We start by summarizing the
continuous-time results. Without loss of generality we assume that $f$ is
normalized such that the Lipschitz constant of the gradient is unity.

######  Theorem 2.

Let $(x(t)$, $u(t))$ be a trajectory satisfying (7) (according to Sec. 3) with
$x(0)\in C$. Let $f$ be $1$-smooth, let $g$ satisfy the Mangasarian-Fromovitz
constraint qualification, and let either $f$ be convex or $2\delta-\beta>0$.
Then, $x(t)$ converges to the set of stationary points, while $u(t)$ converges
to zero. Moreover, each isolated local minimum corresponds to an
asymptotically stable equilibrium in the sense of Lyapunov.

Table 1: The table summarizes convergence rates that arise from different choices of $\alpha$, $\beta$, and $\delta$. Without loss of generality $f$ is normalized, such that the Lipschitz constant of the gradient is unity. variant | $\alpha$ | $\delta$ | $\beta$ | rate $\rho$
---|---|---|---|---
heavy ball | $\sqrt{\mu}$ | $\sqrt{\mu}$ | $0$ | $e^{-\sqrt{\mu}t}$
Nesterov constant parameters | $\sqrt{\mu}-\mu/2$ | $\frac{\sqrt{\mu}}{1+\sqrt{\mu}}$ | $\frac{1-\sqrt{\mu}}{1+\sqrt{\mu}}$ | $e^{-(\sqrt{\mu}-\mu/2)t}$
Nesterov varying parameters | $\frac{2}{t+3}$ | $\frac{3}{2(t+3)}$ | $\frac{t}{t+3}$ | $\frac{9}{(t+3)^{2}}$

The proof of Thm. 2 is included in App. D.

The following proposition demonstrates that the use of momentum combined with
well-chosen damping parameters indeed results in accelerated convergence rates
($\mathcal{O}(1/t^{2})$ in the smooth and convex case, and $e^{-\sqrt{\mu}t}$
in the smooth and strongly convex case).555Continuous-time rates are indeed
meaningful in this context, since both the Lipschitz constant of $\nabla f$
and the constant in front of $\nabla f$ in (7) are fixed to unity. This fixes
the time scale, as any reparametrization of time, i.e., $t=\tau(s)$, where
$\tau:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is a diffeomorphism,
would alter the way $\nabla f$ enters (7). We refer the reader to Muehlebach
and Jordan [8] for a more general discussion.

######  Theorem 3.

Let $C$ be convex and $f$ be $1$-smooth and either convex or strongly convex
with strong convexity constant $\mu>0$. Let the parameters $\alpha$, $\beta$,
$\delta$, and $\rho$ be chosen according to Table 1 and assume that Slater’s
condition holds. Then, for any $x(0)\in\mathbb{R}^{n}$, $u(0)=0$, the
following holds:

$\min\\{0,g(x(0))\\}^{\mathsf{\tiny{T}}}\lambda^{*}\rho(t)\leq\\\
f(x(t))-f(x^{*})\leq\left(\frac{\alpha^{2}}{2}|x(0)-x^{*}|^{2}+f(x(0))-f(x^{*})\right)\rho(t),$

where $x^{*}$ is the minimizer of (1) and $\lambda^{*}$ is a multiplier that
satisfies the Karush-Kuhn-Tucker conditions.

We also demonstrate convergence of the discrete algorithm (8) in a nonconvex
and possibly stochastic setting. For simplicity we state and prove the
deterministic result when $\epsilon=0$ (as becomes apparent from the proof,
the stochastic case with bounded zero-mean gradient perturbations follows from
the same arguments).

######  Theorem 4.

Let $T_{k}=T_{0}/k^{s}$, $k=1,2,\dots$, for some $T_{0}>0$ and $s\in(1/2,1)$,
and let the function $\min\\{0,g_{1}(x)\\}$ have compact level sets. Let $f$
be 1-smooth and either convex or such that $2\delta-\beta>0$, let
$x_{k},u_{k}$ be the iterates defined in (8) with $\epsilon=0$ and the initial
values $(x_{0},u_{0})\in\mathbb{R}^{2n}$ arbitrary, and let $g$ satisfy the
Mangasarian-Fromovitz constraint qualification. If $u_{k}$ is bounded and $f$
has isolated stationary points, then $x_{k}$ converges to a stationary point
of (1), while $u_{k}$ converges to zero.

We note that the restriction $1/2<s<1$ can be loosened to $1/2<s\leq 1$ if
additional assumptions on the damping parameters $\delta$ and $\beta$ are
satisfied (this requires a slightly more detailed proof).

We note that the behavior of algorithm (8) is complex, as it relies on a
_local_ approximation of the feasible set, whereby multiple constraints can
become active or inactive over the course of the optimization. Establishing
Thm. 4 is therefore nontrivial (see the proof in the appendix) and requires
blending ideas from numerical analysis, optimization, and dynamical systems.

We conclude the section by proving an accelerated rate for a slightly modified
version of (8) in a convex setting. The algorithm we will consider is the
following:

$\displaystyle
u_{k+1}=\mathop{\mathrm{argmin}}_{v\in\mathbb{R}^{n}}\frac{1}{2}|v-u_{k}+2\delta
T_{k}u_{k}+\nabla f(y_{k})T_{k}|^{2},$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \text{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
\nabla g_{i}(y_{k})^{\mathsf{\tiny{T}}}v\geq-\alpha
g_{i}(x_{k})-\frac{g_{i}(y_{k})-g_{i}(x_{k})-\nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}\beta u_{k}}{T_{k}},\leavevmode\nobreak\
\leavevmode\nobreak\ i\in[n_{\text{g}}],$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ y_{k}=x_{k}+\beta u_{k},$ $\displaystyle
x_{k+1}=x_{k}+T_{k}u_{k+1}.$ (12)

######  Theorem 5.

Let $f$ be $\mu$-strongly convex and let $g$ be smooth and concave. Let
$l(x)=f(x)-{\lambda^{*}}^{\mathsf{\tiny{T}}}g(x)$ denote the Lagrangian, where
$\lambda^{*}$ denotes an optimal multiplier of (1). Then, the iterates of (12)
with $T_{k}=1/\sqrt{L_{l}}$, $\delta=\sqrt{L_{l}}/(\sqrt{\kappa_{l}}+1)$,
$\beta=T_{k}(1-2\delta T_{k})$, $\alpha=\delta$, and $u_{0}=0$ satisfy:

$\frac{\mu}{2}|x_{k}-x^{*}|^{2}\leq
l(x_{k})-l(x^{*})\leq\left(1-\frac{1}{1+\sqrt{\kappa_{l}}}\right)^{k}\left(\frac{L_{l}}{8}|x_{0}-x^{*}|^{2}+l(x_{0})-l(x^{*})\right),$

for all $k\geq 0$, where $L_{l}$ denotes the smoothness constant of $l$ and
$\kappa_{l}=L_{l}/\mu$.

There are a few differences between (8) and (12). Instead of evaluating
$\nabla g$ at $x_{k}$, the gradient is evaluated at $x_{k}+\beta u_{k}$, which
will simplify our proof (both $\nabla f$ and $\nabla g$ are evaluated at the
same point). Furthermore the right-hand side of the inequality constraint
includes an additional term that can be related to the curvature of $g_{i}$ as
follows

$-\frac{g_{i}(y_{k})-g_{i}(x_{k})-\nabla g_{i}(y_{k})^{\mathsf{\tiny{T}}}\beta
u_{k}}{T_{k}}=\frac{\nabla^{2}g_{i}(\xi_{k})\beta^{2}|u_{k}|^{2}}{2T_{k}},$

for some $\xi_{k}\in\mathbb{R}^{n}$, where Taylor’s theorem has been used in
the second step.

Thm. 5 requires knowledge of the smoothness constant of $l$, which can be
restrictive for certain applications. However, by following the same argument
as in App. A we obtain the following corollary that proves a
$\mathcal{O}(\log(1/\epsilon)/\sqrt{\epsilon})$ convergence rate if $L_{l}$ is
unknown.

######  Corollary 6.

Let the function $f$ be $\mu$-strongly convex and $L$-smooth, let $g$ be
concave and $L_{\text{g}}$-smooth, and let

$B:=\max_{x\in C}|\nabla f(x)|^{2}/(2\mu).$

Then, the iterate $x_{N}$ of (12) with $T_{k}=1/\sqrt{L_{l}}$,
$\delta=\sqrt{L_{l}}/(\sqrt{\kappa_{l}}+1)$, $\beta=T_{k}(1-2\delta T_{k})$,
$\alpha=\delta$, and $L_{l}=L+BL_{\text{g}}/\epsilon$ satisfy
$|x_{N}-x^{*}|\leq\epsilon$, where

$N\\!\geq\\!2\sqrt{\frac{L\\!+\\!BL_{\text{g}}/\epsilon}{\mu}}\\!\left(\\!\\!2\log(1/\epsilon)\\!+\\!\log\left(\frac{5L\\!+\\!5BL_{\text{g}}/\epsilon}{4\mu}\right)\\!+\\!2\log(|x^{*}-x_{0}|)\\!\right)\\!=\\!\mathcal{O}\\!\left(\\!\frac{\log(1/\epsilon)}{\sqrt{\epsilon}}\\!\right),$

and $x^{*}$ denotes the minimizer of
$\min_{x\in\mathbb{R}^{n}}f(x)\leavevmode\nobreak\
\text{s.t.}\leavevmode\nobreak\ g(x)\geq-\epsilon$.

Figure 2: The first panel shows trajectories resulting from (7) (with
parameters $\alpha=0.5,\delta=0.1,\beta=0,\epsilon=0$). The boundaries of
$\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ are highlighted in red. The second
panel shows the results from the discretization (8) with $T_{k}=T=0.1$, while
the third panel shows the results from the discretization (12) with
$T_{k}=T=0.1$. An important difference between (8) and (12) lies in the fact
that only violated constraints are considered in (8), whereas (12) includes
all constraints. This is indicated by the red lines, which denote
$\mathcal{R}_{1}$, $\mathcal{R}_{2}$ in the second panel and
$\gamma_{1}(x,u)=0$, $\gamma_{2}(x,u)=0$ in the third panel.

$0$$50$$100$$150$$200$$250$$300$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$iterations$|f(x_{k})-f^{*}|$Alg.
3 p=1Alg. 4 p=1proj. grad.acc. proj. grad.$1/k^{2}$

Figure 3:  The left panel shows the solution vector of the compressed sensing
problem with $\ell^{1}$ and $\ell^{0.8}$ regularization. The right panel shows
the evolution of the objective function for the different methods. We note
that Alg. 3, Alg. 4, and accelerated projected gradients converge at a similar
rate, which is much faster than gradient descent. We applied the following
settings for Alg. 3 and Alg. 4: $\alpha_{k}=2/(k+3)$, $\delta_{k}=3/(2(k+3))$,
$\beta_{k}=T(1-2\delta_{k}T)$ (see Tab. 1) with $T=1.8$ and $T=2$,
respectively. Accelerated gradient descent corresponds to the algorithm from
Nesterov [36, p. 78, Constant Step Scheme I]. The corresponding trajectories
for Alg. 3 and Alg. 4 for $p<1$ are similar to $p=1$ and are shown in Fig. 4.

$0$$100$$200$$300$$400$$500$$10^{-4}$$10^{-2}$$10^{0}$$10^{2}$iterations$|f(x_{k})-f^{*}|$Alg.
3 $p=0.8$Alg. 4 $p=0.8$$100/k^{2}$

$0$$100$$200$$300$$400$$500$$10^{-3}$$10^{-1}$$10^{1}$iterations$|\nu-\sum_{i=1}^{n}(\bar{x}_{ki})_{\Delta}^{p}|$Alg.
3 $p=0.8$Alg. 4 $p=0.8$$100/k^{2}$

Figure 4:  The figure shows the trajectories of Alg. 3 and Alg. 4 applied to
the compressed sensing problem with $\ell^{0.8}$ regularization. The left
panel shows the evolution of the objective function for the different methods,
whereas the right panel shows the value of the constraint violation. We
applied the following settings for Alg. 3 and Alg. 4: $\alpha_{k}=2/(k+3)$,
$\delta_{k}=3/(2(k+3))$, $\beta_{k}=T(1-2\delta_{k}T)$ (see Tab. 1) with $T=1$
and $\Delta=1e-3$.

##  5 Numerical Examples

The following section is divided into two parts. The first part illustrates
the dynamics of (7) and the discretization via (8) and (12) on a one-
dimensional example and is intended to provide insights concerning the non-
smooth dynamics, as well as the discretization. The second part applies (8)
and (12) to (nonconvex) compressed sensing and large-scale sparse regression
problems. As we will see, our algorithm recovers state-of-the-art performance
for convex relaxations, while also handling nonconvex sparsity constraints in
a seamless manner (traditional projection-based methods cannot be easily
extended to this setting).

###  5.1 Illustrative example

In order to plot trajectories in the phase space we choose $f(x)=(x+2)^{2}/2$
and $g(x)=(x,-x+2)$, where $x$ is scalar. Each constraint $g_{i}(x)\geq 0$ and
its corresponding velocity constraint $\gamma_{i}(x)\geq 0$ induces a region,

$\mathcal{R}_{i}:=\\{(x,u)\in\mathbb{R}^{2}\leavevmode\nobreak\
|\leavevmode\nobreak\ g_{i}(x)\leq 0,\gamma_{i}(x,u)\leq 0\\},$

in the phase space, $i=1,2$, where trajectories are either non-smooth or slide
along the boundary of $\mathcal{R}_{i}$. Outside of $\mathcal{R}_{i}$, the
trajectories follow the smooth motion (6). The first panel in Fig. 2 shows the
continuous-time trajectories (7) along with $\mathcal{R}_{1}$ and
$\mathcal{R}_{2}$. For a given $(x(t_{0}),u(t_{0})^{-})$ an impact happens if
$g_{i}(x(t_{0}))\leq 0$ and $\gamma_{i}(x(t_{0}),u(t_{0}))^{-}<0$, which
ensures that
$\gamma_{i}(x(t_{0}),u(t_{0}))^{+}\geq-\epsilon\gamma_{i}(x(t_{0}),u(t_{0}))^{-}$.
In our example only the case $\gamma_{i}^{+}=-\epsilon\gamma_{i}^{-}$ occurs,
as there are no impacts where more than one constraint participates
($\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ are disjoint). The coefficient of
restitution $\epsilon$ therefore determines the velocity after impact. For
$\epsilon=0$ trajectories end up at the boundary of the set $\mathcal{R}_{i}$,
whereas for $\epsilon>0$ they will leave $\mathcal{R}_{i}$ (in case of
impact). If $g_{i}(x(t_{0}))\leq 0$, $\gamma_{i}(x(t_{0}),u(t_{0}))^{-}=0$, no
impact occurs, ($u(t_{0})=u(t_{0})^{-}=u(t_{0})^{+}$), and trajectories either
leave $\mathcal{R}_{i}$ or slide along its boundary. This depends on the
contribution of the unconstrained dynamics, that is, on the vector
$v_{\text{uc}}(t_{0}):=(u(t_{0}),-2\delta u(t_{0})-\nabla f(x(t_{0})+\beta
u(t_{0}))).$ If $v_{\text{uc}}(t_{0})$ points outwards, trajectories will
leave $\mathcal{R}_{i}$ and follow the unconstrained motion
($\mathrm{d}\lambda_{i}=0$). If $v_{\text{uc}}(t_{0})$ points inwards, there
will be a contribution from
$\mathrm{d}\lambda_{i}=\lambda_{i}(t_{0})\mathrm{d}t$, which ensures that
trajectories slide along the boundary of $\mathcal{R}_{i}$.

The second panel in Fig. 2 shows the trajectories resulting from a
discretization according to (8) with $T_{k}=T=0.1$. We can clearly see the
consequences of including constraints on the velocity level: Trajectories may
become infeasible, since constraints enter (8) only once they are violated.
Nevertheless, even for large time steps $T_{k}=T$ (up to $T\approx 1.8$),
trajectories converge to the unique minimizer of our problem. The third panel
in Fig. 2 shows trajectories from the discretization according to (12) with
$T=0.1$. We note that in (12) each constraint is permanently active (that is,
$I=[n_{\text{g}}]$ in (2)), which explains the contrast between the
trajectories in the second and third panel. In this example the dynamics (12)
are more robust with respect to large time steps, and in fact, convergence can
still be observed for $T=3$.

###  5.2 Nonconvex compressed sensing and image reconstruction

We consider the following $\ell^{p}$-regularized inverse problem:

$\min_{x\in\mathbb{R}^{n}}\frac{1}{2}|Ax-b|^{2}\quad\text{s.t.}\quad|x|_{p}^{p}\leq\nu,$
(13)

where $|x|_{p}$ refers to the $\ell^{p}$ “norm” (we explicitly allow for
$0<p\leq 1$). This has numerous applications in machine learning, statistics,
and signal processing [see, e.g., 37]. The traditional convex approach for
solving such an inverse problem is to set $p=1$ and to leverage the fact that
projections onto the $\ell^{1}$ ball have closed-form solutions. This yields
both accelerated and non-accelerated gradient descent schemes, which in the
setting of (13) are also known under the name of iterative shrinkage-
thresholding (ISTA) and fast iterative shrinkage-thresholding (FISTA) [see,
e.g., 38]. However, when $p<1$, projections onto the $\ell^{p}$ “norm” ball no
longer have closed-form solutions and it is unclear how to generalize
projected gradient algorithms to this setting. Nonetheless the setting $p<1$
can be handled easily with our algorithms, as we highlight next.

In order to handle the absolute value and the $p$th power, we add slack
variables, $\bar{x}\in\mathbb{R}^{n}$ and reformulate (13) as

$\min_{(x,\bar{x})\in\mathbb{R}^{2n}}\frac{1}{2}|Ax-b|^{2}\quad\text{s.t.}\quad-\bar{x}\leq
x\leq\bar{x},\quad\sum_{i=1}^{n}(\bar{x}_{i})_{\Delta}^{p}\leq\nu,$ (14)

where $(\cdot)_{\Delta}^{p}:\mathbb{R}\rightarrow\mathbb{R}$ is continuously
differentiable and approximates $x^{p}$ for $x>0$ ($x^{p}$, $p<1$ is
nondifferentiable at the origin). The approximation $(\cdot)_{\Delta}^{p}$
depends on the approximation parameter $\Delta>0$ and is defined as

$(x)_{\Delta}^{p}:=\begin{cases}x^{p}-\Delta^{p}(1-p)&x\geq\Delta,\\\
p\Delta^{p-1}x,&x<\Delta.\end{cases}$

For $p=1$ our approximation recovers the constraint $|x|_{1}\leq\nu$ exactly
for any $\Delta>0$. The approximation is depicted in Fig. 1 for $p=0.6$ and
$\Delta=0.01$ and highlights that $x^{p}$ is approximated well even for modest
values of $\Delta$. In the numerical examples $\Delta$ is typically set to
$10^{-6}$, which yields an excellent agreement between $x^{p}$ and
$(x)^{p}_{\Delta}$.

Despite the fact that the constraint in (14) is nonlinear and nonconvex for
$p<1$, the optimization in (8) and (12) can be carried out in closed form,
which yields the two algorithms Alg. 3 and Alg. 4 stated in App. C.

#### Nonconvex compressed sensing example:

In the first example, each element of $A\in\mathbb{R}^{100\times 1000}$ is
sampled from a standard normal distribution. The vector $b$ is set to
$Ax^{*}+n/2$, where the components of $n\in\mathbb{R}^{100}$ are sampled from
a standard normal and $x^{*}$ is a vector that contains zeros everywhere
except for 13 randomly chosen entries that are set to one. This gives rise to
a challenging and ill-conditioned optimization problem that includes 1000
decision variables. The left panel in Fig. 3 compares the results computed by
our algorithm for $p=1$ and $p=0.8$, whereas the right panel (solid lines)
compares our approach to projected gradient descent and accelerated projected
gradient descent for $p=1$. We note: i) the quality of the reconstruction for
$p=1$ is significantly worse compared to $p=0.8$ (the parameter $\nu$ was
tuned with five-fold cross validation, which yielded $\nu\approx 13$ in both
cases) and ii) our algorithms decrease the objective function at a similar
rate as accelerated projected gradient for $p=1$. Aside from the computation
of the gradient, all algorithms have the same $\mathcal{O}(n\text{log}(n))$
complexity per iteration, which is determined by a sorting operation that is
used for projecting the iterates onto the $\ell^{1}$ ball or solving the
minimization in (8) or (12), respectively. It is important to note that for
both our algorithms the per-iteration complexity is independent of $p$. Fig. 4
shows the trajectories of Alg. 3 and Alg. 4 for $p=0.8$ and highlights a
convergence rate, both in terms of function value and in terms of constraint
violation, on the order of $1/k^{2}$ also in the nonconvex case.

#### Nonconvex image reconstruction example:

The second example consists of an image reconstruction problem taken from Beck
and Teboulle [38], where $A=RW\in\mathbb{R}^{n\times n},n=65536$, with $R$
representing a Gaussian blur operator, $W$ the inverse of a three-stage Haar
wavelet transform, and $\nu=6\cdot 10^{3}$. The problem given by (14) is of
considerable size and includes 131,072 decision variables and 131,073
constraints. Similar to the previous example, our approach is on par with the
performance of accelerated gradient descent for $p=1$, as shown in Fig. 5.
However, we are also able to solve problems with $p<1$, as highlighted in Fig.
6. Fig. 6 compares the resulting reconstruction of accelerated gradient
descent with $p=1$ compared to our reconstruction $p=0.8$ after 100
iterations, whereby the latter has visibly fewer artifacts. We found that
choice of the damping parameters $\delta_{k}$ and $\beta_{k}$ has a
significant effect on the convergence rate. Our choice
$\delta_{k}=3/(2(k+3)),\beta_{k}=1-2\delta_{k}$, $\alpha_{k}=2/(k+3)$, and
$T=1$ as motivated in Tab. 1, is not optimized and a different schedule might
lead to additional improvements.

Further details, including the derivation of Alg. 3 and Alg. 4 from (8) and
(12), are included in App. C.

$0$$20$$40$$60$$80$$100$$10^{-1}$$10^{1}$$10^{3}$$10^{5}$iterations$|f(x_{k})-f^{*}|$Alg.
3, $p=1$Alg. 4, $p=1$acc. proj. grad.proj. grad.$10^{3}/k^{2}$

$0$$20$$40$$60$$80$$100$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$iterations$|\nu-\sum_{i=1}^{n}(\bar{x}_{ki})_{\Delta}^{p}|$Alg.
3, $p=1$Alg. 4, $p=1$$10^{4}/k^{2}$

Figure 5:  The figure on the left shows the decrease in the objective function
as a function of the iterations for the different algorithms. We note that
Alg. 3, Alg. 4, and accelerated projected gradients converge at a similar
rate, which is substantially faster than gradient descent. We applied the
following settings for Alg. 3 and Alg. 4: $\alpha_{k}=2/(k+3)$,
$\delta_{k}=3/(2(k+3))$, $\beta_{k}=T(1-2\delta_{k}T)$ (see Tab. 1) with
$T=1$. Accelerated gradient descent corresponds to the algorithm by Nesterov
[36, p. 78, Constant Step Scheme I]. The figure on the right shows how
constraint violations decrease as a function of the number of iterations. The
black dashed line indicates a rate of $\mathcal{O}(1/k^{2})$ as a reference.
The corresponding trajectories of Alg. 3 and Alg. 4 for $p<1$ are similar to
$p=1$.

Figure 6: From left to right: Original image, blurred and noisy image, output
after 100 iterations with accelerated gradient descent ($p=1$), output after
100 iterations with Alg. 4 ($p=0.8$). Both the results with $p=1$ and $p=0.8$
lead to reasonable reconstructions, however, $p=0.8$ has fewer artifacts.

##  6 Conclusion

We have introduced a new type of accelerated optimization algorithm for
constrained optimization problems. By imposing constraints on velocities,
rather than on positions, the algorithms avoid projections or optimizations
over the entire feasible set at each iteration. This not only has the
potential to reduce execution time compared to Frank-Wolfe or projected
gradient schemes, but more importantly, expands the range of potential
applications, as constraints are not necessarily required to be convex or to
have a simple structure. We have highlighted important analogies to non-smooth
dynamical systems, and characterized the algorithm’s behavior in continuous
and discrete time.

Acknowledgments Michael Muehlebach thanks the German Research Foundation and
the Branco Weiss Fellowship, administered by ETH Zurich, for the support. We
also acknowledge support from the European Research Council in the form of an
ERC Synergy Grant.

##  Appendix A Supplementary Material of Sec. 2

The application of Thm. 1 requires knowledge of the smoothness constant
$L_{l}$ of $l$, which may require a bound on the optimal dual multiplier
$\lambda^{*}$ if $g$ is nonlinear. This can be viewed as a quantitative
constraint qualification, since the multiplier $\lambda^{*}$ can be bounded by

$|\lambda^{*}|\leq\sup_{x\in C}\frac{|W(x)^{\mathsf{\tiny{T}}}\nabla
f(x)|}{\underline{\sigma}(W(x))},\quad W(x):=(\nabla g_{i}(x))_{i\in I_{x}},$

where $\underline{\sigma}$ denotes the minimum singular value and $W(x)$ the
matrix with columns $\nabla g_{i}(x)$, $i\in I_{x}$.

If such a quantitative constraint qualification is not available, the
following corollary of Thm. 1 still applies and characterizes the convergence
of Alg. 1:

######  Corollary 7.

Let the function $f$ be $\mu$-strongly convex and $L$-smooth, let $g$ be
concave and $L_{\text{g}}$-smooth, and let

$B:=\max_{x\in C}|\nabla f(x)|^{2}/(2\mu).$

Then, the iterate $x_{N}$ of Alg. 1 with $T=1/(L+BL_{g}/\epsilon)$,
$\alpha=\mu$ satisfies $|x_{N}-x^{*}|\leq\epsilon$, where

$N\geq\frac{L+BL_{\text{g}}/\epsilon}{\mu}\left(2\log(1/\epsilon)+\log(\frac{L+BL_{\text{g}}/\epsilon}{\mu})+2\log(|x^{*}-x_{0}|)\right)=\mathcal{O}\left(\frac{\log(1/\epsilon)}{\epsilon}\right)$

and $x^{*}$ denotes the minimizer of

$\min_{x\in\mathbb{R}^{n}}f(x)\quad\text{s.t.}\quad g(x)\geq-\epsilon.$

###### Proof.

We start by proving the following claim:

Claim: Let $\lambda^{*}$ be an optimal multiplier of
$\min_{x\in\mathbb{R}^{n}}f(x)\leavevmode\nobreak\
\text{s.t.}\leavevmode\nobreak\ g(x)\geq-\epsilon$. Then $\lambda^{*}$
satisfies $|\lambda^{*}|_{1}\leq B/\epsilon$.

Proof of the claim: We consider the slightly modified problem

$\min_{x\in\mathbb{R}^{n},\xi\geq
0}f(x)+\bar{B}\xi/\epsilon\quad\text{s.t.}\quad g(x)\geq-\xi,$ (15)

and denote its unique minimizer by $x^{*}(\bar{B}),\xi^{*}(\bar{B})$ (strong
convexity concludes that the minimizer is unique). In the above minimization
the variable $\xi$ can be interpreted as slackness, whereby the factor
$\bar{B}/\epsilon$ penalizes the magnitude of $\xi$. The larger the factor
$\bar{B}/\epsilon$, the smaller the amount of constraint violation by the
resulting $x^{*}(\bar{B})$. The modified problem (15) is motivated by the fact
that corresponding Lagrange dual problem is given by

$\max_{\lambda\geq 0,|\lambda|_{1}\leq\bar{B}/\epsilon}\leavevmode\nobreak\
\min_{x\in\mathbb{R}^{n}}f(x)-\lambda^{\mathsf{\tiny{T}}}g(x),$

which implies that the minimizer $x^{*}(\bar{B})$ satisfies

$x^{*}(\bar{B})=\mathop{\mathrm{argmin}}_{x\in\mathbb{R}^{n}}f(x)-\lambda^{*}(\bar{B})^{\mathsf{\tiny{T}}}g(x),\quad\text{with}\quad|\lambda^{*}(\bar{B})|_{1}\leq\bar{B}/\epsilon,$
(16)

and where $\lambda^{*}(\bar{B})$ is an optimal dual multiplier. The maximum
theorem implies that $x^{*}(\bar{B})$ and $\xi^{*}(\bar{B})$ are continuous
functions of $\bar{B}$. We will show next that for $\bar{B}=B$ the amount of
constraint violation $g(x^{*}(B))$ is guaranteed to be below $\epsilon$.
Combined with the continuity of $\xi^{*}(\bar{B})$, this implies that there
exists some $B^{\prime}\leq B$ for which the corresponding $x^{*}(B^{\prime})$
satisfies
$x^{*}(B^{\prime})=\mathop{\mathrm{argmin}}_{x\in\mathbb{R}^{n},g(x)\geq-\epsilon}f(x)$.
This in turn proves the claim in view of (16).

In order to bound the amount of constraint violation for $\bar{B}=B$, we note
that the following holds for any $x\in C$

$\displaystyle f(x^{*}(\bar{B}))+\bar{B}\xi^{*}(\bar{B})/\epsilon\leq f(x)$
$\displaystyle\leq f(\hat{x})+|\nabla f(x)|^{2}/(2\mu)$ $\displaystyle\leq
f(x^{*}(\bar{B}))+|\nabla f(x)|^{2}/(2\mu),$

where strong convexity of $f$ has been used in the second step and $\hat{x}$
denotes the unconstrained minimizer of $f$. Maximizing the right-hand side
over $x\in C$ and rearranging terms yields

$\xi^{*}(\bar{B})\leq\epsilon B/\bar{B},$

which establishes that $\xi^{*}(B)\leq\epsilon$ and proves the claim.

We now turn to the proof of the corollary. As a result of the claim we
conclude that the minimizer $x^{*}$ satisfies
$x^{*}=\mathop{\mathrm{argmin}}_{x\in\mathbb{R}^{n}}f(x)-{\lambda^{*}}^{\mathsf{\tiny{T}}}g(x)$,
where the $\ell^{1}$ norm of $\lambda^{*}$ is bounded by $B/\epsilon$. As a
result of the choice of $T$ and the parameter $\alpha$, we conclude from Thm.
1 that

$\frac{\mu}{2}|x_{N}-x^{*}|^{2}\leq\left(1-\frac{L+BL_{\text{g}}/\epsilon}{\mu}\right)^{N}\leavevmode\nobreak\
\leavevmode\nobreak\ \frac{L+BL_{\text{g}}/\epsilon}{2}\leavevmode\nobreak\
\leavevmode\nobreak\ |x_{0}-x^{*}|^{2},$

where the smoothness of $f(x)-{\lambda^{*}}^{\mathsf{\tiny{T}}}g(x)$ has been
used to upper bound the right-hand side. The corollary follows from using the
identity
$1-(L+BL_{\text{g}}/\epsilon)/\mu\leq\exp(-(L+BL_{\text{g}}/\epsilon)/\mu)$,
which yields

$|x_{N}-x^{*}|^{2}\leq\frac{L+BL_{\text{g}}/\epsilon}{\mu}|x_{0}-x^{*}|^{2}\\\
\exp\left(-2\log(1/\epsilon)-\log(\frac{L+BL_{\text{g}}/\epsilon}{\mu})-2\log(|x^{*}-x_{0}|)\right)=\epsilon^{2}$

and proves the desired result. ∎

The section concludes with the proof of Thm. 1.

###### Proof.

The proof hinges on the fact that $\alpha(x^{*}-x_{k})\in V_{\alpha}(x_{k})$,
which follows from concavity of $g$. We define $v_{k}:=(x_{k+1}-x_{k})/T$ to
be the velocity. From the fact that $v_{k}=\mathop{\mathrm{argmin}}_{v\in
V_{\alpha}(x_{k})}|v+\nabla f(x_{k})|^{2}$ we conclude

$\displaystyle|v_{k}+\nabla f(x_{k})|^{2}/2\leq|\alpha(x^{*}-x_{k})+\nabla
f(x_{k})|^{2}/2,$

which can be rearranged to

$\displaystyle|v_{k}|^{2}/2+\nabla f(x_{k})^{\mathsf{\tiny{T}}}v_{k}$
$\displaystyle\leq\alpha^{2}|x^{*}-x_{k}|^{2}/2+\alpha\nabla
f(x_{k})^{\mathsf{\tiny{T}}}(x^{*}-x_{k})$
$\displaystyle\leq(\alpha^{2}-\alpha\mu)|x^{*}-x_{k}|^{2}/2+\alpha(f(x^{*})-f(x_{k}))$
$\displaystyle\leq\alpha(f(x^{*})-f(x_{k})),$ (17)

where we have used the strong convexity of $f$ in the second step and the fact
that $\alpha\leq\mu$ in the third step. Next we use the smoothness of $l$,
which yields

$\displaystyle l(x_{k+1})-l(x_{k})$ $\displaystyle\leq T\nabla
l(x_{k})^{\mathsf{\tiny{T}}}v_{k}+\frac{T^{2}L_{l}}{2}|v_{k}|^{2}$
$\displaystyle\leq T\nabla
f(x_{k})^{\mathsf{\tiny{T}}}v_{k}-T{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g(x_{k})^{\mathsf{\tiny{T}}}v_{k}+\frac{T^{2}L_{l}}{2}|v_{k}|^{2}$
$\displaystyle\leq\alpha
T(f(x^{*})-f(x_{k}))-T{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g(x_{k})^{\mathsf{\tiny{T}}}v_{k}+T\frac{TL_{l}-1}{2}|v_{k}|^{2},$

where we used inequality (17) in the third step. We further note that by
definition of $V_{\alpha}(x_{k})$, $v_{k}\in V_{\alpha}(x_{k})$ satisfies
$\nabla g(x_{k})^{\mathsf{\tiny{T}}}v_{k}\geq-\alpha g(x_{k})$ and therefore
$-{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g(x_{k})^{\mathsf{\tiny{T}}}v_{k}\leq-\alpha{\lambda^{*}}^{\mathsf{\tiny{T}}}g(x_{k})$.
In addition, complementary slackness implies
${\lambda^{*}}^{\mathsf{\tiny{T}}}g(x^{*})=0$ and therefore

$l(x_{k+1})-l(x_{k})\leq\alpha
T(l(x^{*})-l(x_{k}))+T\frac{TL_{l}-1}{2}|v_{k}|^{2}.$

Subtracting and adding $l(x^{*})$ on the left-hand side and using $T\leq
1/L_{l}$ yields

$\frac{\mu}{2}|x_{k}-x^{*}|^{2}\leavevmode\nobreak\ \leavevmode\nobreak\
\leq\leavevmode\nobreak\ \leavevmode\nobreak\
l(x_{k})-l(x^{*})\leavevmode\nobreak\ \leavevmode\nobreak\
\leq\leavevmode\nobreak\ \leavevmode\nobreak\ (1-\alpha
T)^{k}\leavevmode\nobreak\ \leavevmode\nobreak\ (l(x_{0})-l(x^{*})),$

which concludes the proof. ∎

##  Appendix B Pseudo-code of algorithm (8)

This section lists the pseudo-code of the resulting discrete-time algorithm
(8).

parameters $\beta\geq 0$, $\delta>0$ $\triangleright$ typically below 1; can
be time varying

step size $T_{k}$ $\triangleright$ $1/T_{k}\sim$ smoothness of objective

constants $\alpha>0$, $0\leq\epsilon\leq 1$ $\triangleright$ $\alpha$ such
that $\alpha T_{k}\leq 1$

initial condition $x_{0}\in\mathbb{R}^{n}$, $u_{0}\in\mathbb{R}^{n}$
$\triangleright$ usually $u_{0}=0$

$\epsilon_{\text{const}}\geq 0$ $\triangleright$ constraint violation
tolerance

$x_{k}\leftarrow x_{0}$, $u_{k}\leftarrow u_{0}$

for $k=0,1,\dots$ do

$r_{k}\leftarrow u_{k}-2\delta T_{k}u_{k}-\nabla f(x_{k}+\beta u_{k})T_{k}$
$\triangleright$ unconstrained update

$w\leftarrow\\{\\}$, $W\leftarrow\\{\\}$

for $i=1,\dots,n_{\text{g}}$ do $\triangleright$ check constraint violations

if $g_{i}(x_{k})\leq\epsilon_{\text{const}}$ then

$w\leftarrow(w,-\alpha g_{i}(x_{k})-\epsilon\min\\{\nabla
g_{i}(x_{k})^{\mathsf{\tiny{T}}}u_{k}+\alpha g_{i}(x_{k}),0\\})$

$W\leftarrow(W,\nabla g_{i}(x_{k})^{\mathsf{\tiny{T}}})$ $\triangleright$ add
violated constraint

end if

end for

$u_{k+1}\leftarrow\mathop{\mathrm{argmin}}_{v\in\mathbb{R}^{n}}|v-r_{k}|^{2}$
s.t. $Wv\geq w$ $\triangleright$ solve (8)

$x_{k}\leftarrow x_{k}+T_{k}u_{k+1}$ $\triangleright$ update positions

$u_{k}\leftarrow u_{k+1}$ $\triangleright$ update velocities

end for

##  Appendix C Additional Details on the Numerical Experiments

We apply (8) and (12) to (14) and start by defining the constraint function
$g$ as $g(x,\bar{x}):=(g_{1}(x,\bar{x}),g_{2}(x,\bar{x}))$, where
$g_{1}(x,\bar{x}):=\bar{x}+x$ and $g_{2}(x,\bar{x}):=\bar{x}-x$. The decision
variables $(x,\bar{x})$ will be represented with the iterates
$x_{k}\in\mathbb{R}^{n}$, $\bar{x}_{k}\in\mathbb{R}^{n}$, and the
corresponding velocities will be denoted by $u_{k}\in\mathbb{R}^{n}$,
$\bar{u}_{k}\in\mathbb{R}^{n}$, where $k$ refers to the iteration number. As a
result, the optimization in (8), which is used to determine the velocities
$(u_{k+1},\bar{u}_{k+1})$ can be expressed as

$\displaystyle(u_{k+1},\bar{u}_{k+1})=$
$\displaystyle\mathop{\mathrm{argmin}}_{(v,\bar{v})\in\mathbb{R}^{2n}}\frac{1}{2}|v-r_{k}|^{2}+\frac{1}{2}|\bar{v}-\bar{r}_{k}|^{2}\quad\text{s.t.}$
(18) $\displaystyle\quad\bar{v}_{i}+v_{i}+\tilde{g}_{1i}\geq 0,\quad\forall
i\in I_{1k},\quad\bar{v}_{j}-v_{j}+\tilde{g}_{2j}\geq 0,\quad\forall j\in
I_{2k},$ $\displaystyle\quad-\sum_{i=1}^{n}w_{l}\bar{v}_{l}+\tilde{g}_{3}\geq
0,$

where the last constraint only enters if
$\sum_{i=1}^{n}(\bar{x}_{ki})_{\Delta}^{p}\geq\nu$ holds, and where $I_{ik}$
denotes the set of constraints in $g_{i}$ that are active at iteration $k$,
that is,

$I_{1k}:=\\{i\in\mathbb{Z}\leavevmode\nobreak\ |\leavevmode\nobreak\
g_{1i}(x_{k},\bar{x}_{k})\leq 0\\},\quad
I_{2k}:=\\{i\in\mathbb{Z}\leavevmode\nobreak\ |\leavevmode\nobreak\
g_{2i}(x_{k},\bar{x}_{k})\leq 0\\}.$

The variables $\tilde{g}_{1}$, $\tilde{g}_{2}$, $\tilde{g}_{3}$, and $w$ are
defined as follows: $\tilde{g}_{1}:=\alpha g_{1}(x_{k},\bar{x}_{k})$,
$\tilde{g}_{2}:=\alpha g_{2}(x_{k},\bar{x}_{k})$,
$\tilde{g}_{3}:=\alpha(\nu-\sum_{i=1}^{n}(\bar{x}_{ki})_{\Delta}^{p})$, and
$w_{i}:=\text{d}(x)^{p}_{\Delta}/\text{d}x|_{x=\bar{x}_{ki}}$. We have further
set the constant $\epsilon=0$ when deriving (18) in order to simplify the
exposition. Moreover, $r_{k}\in\mathbb{R}^{n}$ and
$\bar{r}_{k}\in\mathbb{R}^{n}$ are given by

$r_{k}:=u_{k}-2\delta Tu_{k}-A^{\mathsf{\tiny{T}}}(A(x_{k}+\beta
u_{k})-b)T,\qquad\bar{r}_{k}:=\bar{u}_{k}-2\delta T\bar{u}_{k},$

where $\delta>0$ and $\beta\geq 0$ are damping parameters and $T>0$ denotes
the step size. In fact, the optimization in (12) has the same form as (18)
with slightly different definitions of $\tilde{g}_{1},\tilde{g}_{2},$ and
$\tilde{g}_{3}$ and all constraints being permanently active, which means, in
particular, $I_{1k}=\\{1,\dots,n\\}$ and $I_{2k}=\\{1,\dots,n\\}$. In both
cases, we apply the following change of variables
$(v,\bar{v})\rightarrow(\xi,\bar{\xi})$,

$v=\xi-\bar{\xi}-\frac{\tilde{g}_{1}-\tilde{g}_{2}}{2},\qquad\bar{v}=\xi+\bar{\xi}-\frac{\tilde{g}_{1}+\tilde{g}_{2}}{2},$

which transforms (18) to a projection onto a weighted simplex. The projection
onto a weighted simplex has a closed-form solution and takes the form

$\displaystyle\mathop{\mathrm{argmin}}_{(\xi,\bar{\xi})\in\mathbb{R}^{2n}}$
$\displaystyle\frac{1}{2}|\xi-
q_{1}|^{2}+\frac{1}{2}|\bar{\xi}-q_{2}|^{2}\quad\text{s.t.}$ (19)
$\displaystyle\xi_{i}\geq 0,\leavevmode\nobreak\ \forall i\in
I_{1k},\quad\bar{\xi}_{j}\geq 0,\leavevmode\nobreak\ \forall j\in
I_{2k},\quad\sum_{i=1}^{n}w_{i}\xi_{i}+\sum_{i=1}^{n}w_{i}\bar{\xi}_{i}\leq\bar{\nu},$

where the last constraint only enters if
$\sum_{i=1}^{n}(\bar{x}_{ki})_{\Delta}^{p}\geq\nu$ and where $q_{1}$, $q_{2}$,
and $\bar{\nu}$ arise from the change of variables. The weights $w_{i}$ are
given by $\text{d}(x)_{\Delta}^{p}/\text{d}x$ and are therefore nonnegative by
construction of $(\cdot)_{\Delta}^{p}$. The closed-form solution of (19) is
summarized in Alg. 2 below. The solution hinges on the fact that the
constraint
$\sum_{i=1}^{n}w_{i}\xi_{i}+\sum_{i=1}^{n}w_{i}\bar{\xi}_{i}\leq\bar{\nu}$ can
be included in the objective function by means of a dual multiplier, at which
point the resulting minimization corresponds to a simple projection on
$\mathbb{R}_{\geq 0}$ for each $i\in I_{1k}$ and $i\in I_{2k}$. The optimal
multiplier is computed by an iterative procedure that is based on a sorting
operation. The worst-case complexity of Alg. 2 is therefore
$\mathcal{O}(n\text{log}(n))$.

The resulting algorithms for solving the minimization in (8) and (12) are
summarized in Alg. 3 and Alg. 4, respectively.

Algorithm 2 Projection onto the weighted simplex:
$\mathop{\mathrm{argmin}}_{\xi\in\mathbb{R}^{n}}|\xi-q|^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
\xi_{i}\geq 0,\forall i\in I,\leavevmode\nobreak\ \leavevmode\nobreak\
w^{\mathsf{\tiny{T}}}\xi\leq\bar{\nu}$

$q\in\mathbb{R}^{n}$, index set $I$, weight $w\in\mathbb{R}^{n}$, and bound
$\bar{\nu}\in\mathbb{R}$

$w\geq 0$ and $\bar{\nu}\geq 0$

$\xi_{i}\leftarrow\begin{cases}0,&\text{if}\leavevmode\nobreak\ i\in
I\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ q_{i}<0,\\\
q_{i},&\text{else}.\end{cases}$ $\triangleright$ suppose constraint
$w^{\mathsf{\tiny{T}}}\xi\leq\bar{\nu}$ is not active

if $w^{\mathsf{\tiny{T}}}\xi>\bar{\nu}$ then $\triangleright$ constraint
$w^{\mathsf{\tiny{T}}}\xi\leq\bar{\nu}$ is active

$q_{\text{s}}/w_{\text{s}}\leftarrow\text{sort}(\\{q_{i}/w_{i}\\}_{i\in
I},\text{descending})$ $\triangleright$ sort s.t.
$q_{\text{s}1}/w_{\text{s}1}\geq q_{\text{s}2}/w_{\text{s}2}\geq\dots$

$m\leftarrow 1$$\triangleright$ compute multiplier $\lambda$ for
$w^{\mathsf{\tiny{T}}}\xi\leq\bar{\nu}$

while $m\leq|I|$ do

if $\sum_{i\not\in
I}w_{i}q_{i}+\sum_{j<m}w_{\text{s}j}q_{\text{s}j}-(\sum_{i\not\in
I}w_{i}^{2}+\sum_{j<m}w_{\text{s}j}^{2})q_{\text{s}m}/w_{\text{s}m}>\bar{\nu}$
then

break

end if

$m\leftarrow m+1$

end while

$\lambda\leftarrow(\sum_{i\not\in
I}w_{i}q_{i}+\sum_{j<m}w_{\text{s}j}q_{\text{s}j}-\bar{\nu})/(\sum_{i\not\in
I}w_{i}^{2}+\sum_{j<m}w_{\text{s}j}^{2})$

$\xi\leftarrow q-w\lambda$ $\triangleright$ given $\lambda$, perform proj. on
$\mathbb{R}_{\geq 0}$ for each $i\in I$

$\xi_{i}\leftarrow 0$ if $i\in I$ and $\xi_{i}<0$

end if

return $\xi$

Algorithm 3 Algorithm (8) applied to (13).

damping parameters $\beta_{k}\geq 0$, $\delta_{k}>0$; step size $T_{k}$;
$L=|A|^{2}$; $\alpha_{k}>0$

initial condition $x_{0}\in\mathbb{R}^{n}$; approximation parameter $\Delta>0$

$x_{k}\leftarrow x_{0},\bar{x}_{k}\leftarrow\text{abs}(x_{0}),u_{k}\leftarrow
0,\bar{u}_{k}\leftarrow 0$ $\triangleright$ different initialization is also
possible

for $k=0,1,\dots$ do

$r_{k}\leftarrow
u_{k}-2\delta_{k}T_{k}u_{k}-T_{k}A^{\mathsf{\tiny{T}}}(A(x_{k}+\beta_{k}u_{k})-b)/L$
$\triangleright$ unconstrained update

$\bar{r}_{k}\leftarrow\bar{u}_{k}-2\delta T_{k}\bar{u}_{k}$

$\tilde{g}_{1}\leftarrow\alpha_{k}(x_{k}+\bar{x}_{k})$ $\triangleright$
evaluate constraints

$\tilde{g}_{2}\leftarrow\alpha_{k}(-x_{k}+\bar{x}_{k})$

$\tilde{g}_{3}\leftarrow\alpha_{k}(\nu-\sum_{i=1}^{n}(\bar{x}_{ki})_{\Delta}^{p})$

if $\tilde{g}_{3}>0$ then $\triangleright$ evaluate gradient of $\ell^{p}$
constraint

$w\leftarrow 0$

else

$w_{i}\leftarrow\text{d}(x)^{p}_{\Delta}/\text{d}x|_{x=\bar{x}_{ki}},\quad
i=1,2,\dots,n,$

end if

$q\leftarrow(\tilde{g}_{1}/2+(r_{k}+\bar{r}_{k})/2,\leavevmode\nobreak\
\leavevmode\nobreak\ \tilde{g}_{2}/2+(-r_{k}+\bar{r}_{k})/2)$ $\triangleright$
change of variables

$\bar{\nu}\leftarrow\tilde{g}_{3}+w^{\mathsf{\tiny{T}}}(\tilde{g}_{1}+\tilde{g}_{2})/2$

$I\leftarrow\\{i\in\mathbb{Z}:(\tilde{g}_{1},\tilde{g}_{2})_{i}\leq 0\\}$
$\triangleright$ active constraints

$(\xi,\bar{\xi})\leftarrow\text{ProjectionOnWeightedSimplex}(q,I,(w,w),\bar{\nu})$
$\triangleright$ apply Alg. 2

$u_{k}\leftarrow\xi-\bar{\xi}-(\tilde{g}_{1}-\tilde{g}_{2})/2$$\triangleright$
update velocity

$\bar{u}_{k}\leftarrow\xi+\bar{\xi}-(\tilde{g}_{1}+\tilde{g}_{2})/2$

$x_{k}\leftarrow x_{k}+T_{k}u_{k}$ $\triangleright$ update position

$\bar{x}_{k}\leftarrow\bar{x}_{k}+T_{k}\bar{u}_{k}$

end for

Algorithm 4 Algorithm (12) applied to (13).

damping parameters $\beta_{k}\geq 0$, $\delta_{k}>0$; step size $T_{k}$;
$L=|A|^{2}$; $\alpha_{k}>0$

initial condition $x_{0}\in\mathbb{R}^{n}$; approximation parameter $\Delta>0$

$x_{k}\leftarrow x_{0},\bar{x}_{k}\leftarrow\text{abs}(x_{0}),u_{k}\leftarrow
0,\bar{u}_{k}\leftarrow 0$ $\triangleright$ different initialization is also
possible

for $k=0,1,\dots$ do

$r_{k}\leftarrow
u_{k}-2\delta_{k}T_{k}u_{k}-T_{k}A^{\mathsf{\tiny{T}}}(A(x_{k}+\beta_{k}u_{k})-b)/L$
$\triangleright$ unconstrained update

$\bar{r}_{k}\leftarrow\bar{u}_{k}-2\delta T_{k}\bar{u}_{k}$

$w_{i}\leftarrow\text{d}(x)^{p}_{\Delta}/\text{d}x|_{x=\bar{x}_{ki}+\beta_{k}\bar{u}_{ki}},\quad
i=1,2,\dots,n,$ $\triangleright$ gradient of $\ell^{p}$ constraint

$g_{1}\leftarrow x_{k}+\bar{x}_{k}$ $\triangleright$ evaluate constraints at
$x_{k}$

$g_{2}\leftarrow-x_{k}+\bar{x}_{k}$

$g_{3}\leftarrow\nu-\sum_{i=1}^{n}(\bar{x}_{ki})_{\Delta}^{p}$

$g_{y1}\leftarrow
x_{k}+\beta_{k}u_{k}+\bar{x}_{k}+\beta_{k}\bar{u}_{k}$$\triangleright$
evaluate constraints at $x_{k}+\beta_{k}u_{k}$

$g_{y2}\leftarrow-(x_{k}+\beta_{k}u_{k})+\bar{x}_{k}+\beta_{k}\bar{u}_{k}$

$g_{y3}\leftarrow\nu-\sum_{i=1}^{n}(\bar{x}_{ki}+\beta_{k}\bar{u}_{ki})_{\Delta}^{p}$

$\tilde{g}_{1}\leftarrow\alpha_{k}g_{1}+(g_{y1}-g_{1}-\beta_{k}(u_{k}+\bar{u}_{k}))/T_{k}$
$\triangleright$ r.h.s. of constraints in (12)

$\tilde{g}_{2}\leftarrow\alpha_{k}g_{2}+(g_{y2}-g_{2}-\beta_{k}(-u_{k}+\bar{u}_{k}))/T_{k}$

$\tilde{g}_{3}\leftarrow\alpha_{k}g_{3}+(g_{y3}-g_{3}+\beta_{k}w^{\mathsf{\tiny{T}}}\bar{u}_{k})/T_{k}$

$q\leftarrow(\tilde{g}_{1}/2+(r_{k}+\bar{r}_{k})/2,\leavevmode\nobreak\
\leavevmode\nobreak\ \tilde{g}_{2}/2+(-r_{k}+\bar{r}_{k})/2)$ $\triangleright$
change of variables

$\bar{\nu}\leftarrow\tilde{g}_{3}+w^{\mathsf{\tiny{T}}}(\tilde{g}_{1}+\tilde{g}_{2})/2$

$(\xi,\bar{\xi})\leftarrow\text{ProjectionOnWeightedSimplex}(q,\\{1,\dots,2n\\},(w,w),\bar{\nu})$
$\triangleright$ apply Alg. 2

$u_{k}\leftarrow\xi-\bar{\xi}-(\tilde{g}_{1}-\tilde{g}_{2})/2$$\triangleright$
update velocity

$\bar{u}_{k}\leftarrow\xi+\bar{\xi}-(\tilde{g}_{1}+\tilde{g}_{2})/2$

$x_{k}\leftarrow x_{k}+T_{k}u_{k}$ $\triangleright$ update position

$\bar{x}_{k}\leftarrow\bar{x}_{k}+T_{k}\bar{u}_{k}$

end for

##  Appendix D Proof of Thm. 2

From the analysis of the smooth motion in Sec. 3, we infer that
$\gamma_{i}(x(t),u(t))\geq 0$ for all $t$ such that $i\in I_{x(t)}$ (a.e.). By
Grönwall’s inequality, this implies $g(x(t))\geq 0$ for all $t\geq 0$ and
therefore $x(t)\in C$ for all $t\in[0,\infty)$.

We consider the following Lyapunov function

$V(x,u)=\frac{1}{2}|u|^{2}+f(x),$ (20)

which is bounded below by assumption. We investigate how $V(x(t),u(t))$
evolves along the trajectories of (7). The differential measure corresponding
to $V(x(t),u(t))$ can be expressed as [33]:

$\displaystyle\mathrm{d}V$
$\displaystyle=\frac{1}{2}(u^{+}+u^{-})^{\mathsf{\tiny{T}}}\mathrm{d}u+\nabla
f(x)^{\mathsf{\tiny{T}}}u\leavevmode\nobreak\ \mathrm{d}t$
$\displaystyle=-2\delta|u|^{2}\leavevmode\nobreak\ \mathrm{d}t-(\nabla
f(x+\beta u)-\nabla f(x))^{\mathsf{\tiny{T}}}u\leavevmode\nobreak\
\mathrm{d}t+\sum_{i\in
I_{x}}\left(\frac{1}{2}(\gamma_{i}^{+}+\gamma_{i}^{-})-\alpha
g_{i}\right)\mathrm{d}\lambda_{i}.$

The second line follows from replacing $\mathrm{d}u$ with (7), using the fact
that the Lebesgue measure captures only the smooth motion, which means, for
example, $(u^{+})^{\mathsf{\tiny{T}}}u\leavevmode\nobreak\
\mathrm{d}t=(u^{-})^{\mathsf{\tiny{T}}}u\leavevmode\nobreak\
\mathrm{d}t=|u|^{2}\mathrm{d}t$, and adding and subtracting $\alpha
g_{i}\mathrm{d}\lambda_{i}$. From the assumption $2\delta-\beta>0$ (or
alternatively by convexity of $f$) we can upper bound $\mathrm{d}V$ by

$\mathrm{d}V\leq-c_{2}|u|^{2}\mathrm{d}t+\sum_{i\in
I_{x}}\left(\frac{1}{2}(\gamma_{i}^{+}+\gamma_{i}^{-})-\alpha
g_{i}\right)\mathrm{d}\lambda_{i},$

where $c_{2}>0$ is constant. The summand in the second part of the expression
can be rewritten as

$-\alpha
g_{i}\mathrm{d}\lambda_{i}+\frac{1}{2}(\gamma_{i}^{+}+\epsilon\gamma_{i}^{-})\mathrm{d}\lambda_{i}+\frac{1-\epsilon}{2}(\gamma_{i}^{-})\mathrm{d}\lambda_{i}.$

The fact that $g_{i}(x(t))=0$ for all $i\in I_{x(t)}$ ($x(t)$ remains
feasible) implies that the first term vanishes. The second term vanishes due
to the complementarity condition in (7). The third term is guaranteed to be
non-positive, since, on the one hand $\mathrm{d}\lambda_{i}\geq 0$ (see (7)),
and on the other hand, $\gamma_{i}^{-}\leq 0$ in case of impact, and
$\gamma_{i}^{-}\leavevmode\nobreak\
\mathrm{d}\lambda_{i}=\gamma_{i}\mathrm{d}\lambda_{i}=0$ in case of smooth
motion (see again (7)). We therefore conclude that

$\mathrm{d}V\leq-c_{2}|u|^{2}\mathrm{d}t,$ (21)

which means that $V(x(t),u(t))$ is monotonically decreasing in $t$.

For proving convergence to stationary points, we note that since $u(t)$ is of
locally bounded variation with no singular part (by assumption), it can be
decomposed in an absolutely continuous function and a piecewise constant
function [Lebesgue decomposition, see for example Ch. 3 of 33]. As a result,
$\dot{u}(t)$ (whenever it exists) is uniformly locally integrable [see 39].
Combined with the fact that $u$ is square integrable, which follows from (21)
and the observation that $V$ is bounded below, we conclude that
$u(t)\rightarrow 0$ as $t\rightarrow\infty$ by a variant of Barbalat’s lemma
[see 39]. We recall that $x(t)$ is bounded ($C$ is compact) and consider any
cluster point $\bar{x}$ of $x(t)$, which means that there exists a sequence
$t_{k}\rightarrow\infty$ such that $x(t_{k})\rightarrow\bar{x}$. We pick any
$\tau>0$ and consider

$u(t_{k}+\tau)^{+}-u(t_{k})^{-}=\int_{t_{k}}^{t_{k}+\tau}\mathrm{d}u=\int_{t_{k}}^{t_{k}+\tau}-2\delta
u-\nabla f(x+\beta
u)\mathrm{d}t+\sum_{i=1}^{n_{\text{g}}}\int_{t_{k}}^{t_{k}+\tau}\nabla
g_{i}(x)\mathrm{d}\lambda_{i},$

where we set $\mathrm{d}\lambda_{i}=0$ whenever $i\not\in I_{x(t)}$ to
simplify notation. The previous expression is guaranteed to vanish for
$k\rightarrow\infty$, which by continuity of $\nabla g_{i}$ and $\nabla f$
means that

$-\nabla f(\bar{x})\tau+\sum_{i=1}^{n_{\text{g}}}\nabla
g_{i}(\bar{x})\int_{t_{k}}^{t_{k}+\tau}\mathrm{d}\lambda_{i}\rightarrow 0$
(22)

as $k\rightarrow\infty$. We note that by continuity of $g$, there exits a
constant $k_{0}>0$ and $c_{\text{lg}}>0$ such that
$g_{i}(x(t_{k}))>c_{\text{lg}}$ for all $k>k_{0}$ and all $i\not\in
I_{\bar{x}}$. From the fact that $u(t)\rightarrow 0$ we infer that for large
enough $k$, $g_{i}(x(t))>c_{\text{lg}}/2$ for all $t\in[t_{k},t_{k}+\tau]$ and
$i\not\in I_{\bar{x}}$. We introduce the notation

$\int_{t_{k}}^{t_{k}+\tau}\mathrm{d}\lambda_{i}:=\lambda_{k}^{i}\tau,$

for all $k>0$ and all $i\in\\{1,2,\dots,n_{\text{g}}\\}$ and conclude that
$\lambda_{k}^{i}\geq 0$ and, for all $k$ large enough, $\lambda_{k}^{i}=0$ if
$i\not\in I_{\bar{x}}$.

We argue next that $\lambda_{k}^{i}$ is bounded for all
$i\in\\{1,2,\dots,n_{\text{g}}\\}$. We argue by contradiction and consider the
sequence
$(\lambda_{k}^{1},\dots,\lambda_{k}^{n_{\text{g}}})/|(\lambda_{k}^{1},\dots,\lambda_{k}^{n_{\text{g}}})|$,
which is guaranteed to be bounded, even though
$|(\lambda_{k}^{1},\dots,\lambda_{k}^{n_{\text{g}}})|\rightarrow\infty$. Upon
dividing (22) by $|(\lambda_{k}^{1},\dots,\lambda_{k}^{n_{\text{g}}})|$ we
conclude that

$\sum_{i=1}^{n_{\text{g}}}\nabla g_{i}(\bar{x})\xi_{i}=0,$

where $\xi$ denotes an accumulation point of
$(\lambda_{k}^{1},\dots,\lambda_{k}^{n_{\text{g}}})/|(\lambda_{k}^{1},\dots,\lambda_{k}^{n_{\text{g}}})|$,
which satisfies $\xi_{i}\geq 0$, $\xi_{i}=0$ for all $i\not\in I_{\bar{x}}$,
and $|\xi|=1$. However, the fact that $\xi\neq 0$ contradicts the constraint
qualification, since it precludes the existence of a vector $w$ such that
$w^{\mathsf{\tiny{T}}}\nabla g_{i}(\bar{x})>0$ for all $i\in I_{\bar{x}}$
(which would mean $\sum_{i=1}^{n_{\text{g}}}w^{\mathsf{\tiny{T}}}\nabla
g_{i}(\bar{x})\xi_{i}>0$). This implies that $\lambda_{k}^{i}$ is bounded for
all $i\in\\{1,2,\dots,n_{\text{g}}\\}$.

We take any accumulation point of
$(\lambda_{k}^{1},\dots,\lambda_{k}^{n_{\text{g}}})$, which we denote by
$\bar{\lambda}$. The accumulation point $\bar{\lambda}$ satisfies
$\bar{\lambda}\geq 0$, $\bar{\lambda}_{i}=0$ for all $i\not\in I_{\bar{x}}$
(complementary slackness), and by (22), $-\nabla
f(\bar{x})+\sum_{i=1}^{n_{\text{g}}}\nabla g_{i}(\bar{x})\bar{\lambda}_{i}=0$.
Hence, $\bar{x}$ and $\bar{\lambda}$ satisfy the Karush-Kuhn-Tucker conditions
of (1), and $\bar{x}$ is stationary.

We conclude the section by proving asymptotic stability of isolated local
minima (in the sense of Lyapunov). We accordingly pick any isolated local
minimum $x^{*}$ and note that the function $|u|^{2}/2+f(x)-f(x^{*})$ is
positive definite in a neighborhood of $(x^{*},0)$. We conclude from (21) that
$x^{*}$ is therefore stable in the sense of Lyapunov. We have already shown
attractiveness (see previous paragraphs) and therefore conclude that $x^{*}$
is asymptotically stable in the sense of Lyapunov. ∎

##  Appendix E Proof of Thm. 3

#### i) Heavy ball ($\alpha=\sqrt{\mu}$):

We consider the evolution of the function

$W(x,u)=\frac{1}{2}|\alpha(x-x^{*})+u|^{2}+f(x)-f(x^{*}),$

along the trajectories $(x(t),u(t))$. The corresponding differential measure
$\mathrm{d}W$ is given by [33]:

$(\alpha(x-x^{*})+(u^{+}+u^{-})/2)^{\mathsf{\tiny{T}}}(\alpha
u\leavevmode\nobreak\ \mathrm{d}t+\mathrm{d}u)+\nabla
f(x)^{\mathsf{\tiny{T}}}u\leavevmode\nobreak\ \mathrm{d}t.$

By following the same steps as in the proof of Thm. 2 we obtain

$(u^{+}+u^{-})/2^{\mathsf{\tiny{T}}}\mathrm{d}u+\nabla
f(x)^{\mathsf{\tiny{T}}}u\leavevmode\nobreak\
\mathrm{d}t\leq-2\delta|u|^{2}\leavevmode\nobreak\ \mathrm{d}t+\sum_{i\in
I_{x}}-\alpha g_{i}(x)\mathrm{d}\lambda_{i}.$

However, in contrast to Thm. 2, we allow also for infeasible initial
conditions, and therefore the term $-\alpha g_{i}(x)\mathrm{d}\lambda_{i}$
remains. This yields the following upper bound on $\mathrm{d}W$:

$\mathrm{d}W\leq-(2\delta-\alpha)|u|^{2}\mathrm{d}t-\alpha(x-x^{*})^{\mathsf{\tiny{T}}}\nabla
f(x)\mathrm{d}t+(\alpha^{2}-2\alpha\delta)(x-x^{*})^{\mathsf{\tiny{T}}}u\mathrm{d}t\\\
+\sum_{i\in I_{x}}(\alpha\nabla g_{i}(x)^{\mathsf{\tiny{T}}}(x-x^{*})-\alpha
g_{i}(x))\mathrm{d}\lambda_{i}.$ (23)

The fact that $f$ is strongly convex means that the following holds:

$\displaystyle-(x-x^{*})^{\mathsf{\tiny{T}}}\nabla f(x)\leq$
$\displaystyle-(f(x)-f(x^{*}))-\frac{\mu}{2}|x-x^{*}|^{2}.$

In addition, convexity of $C$ implies that each $g_{i}$ is concave, and
therefore

$\nabla g_{i}(x)^{\mathsf{\tiny{T}}}(x^{*}-x)\geq g_{i}(x^{*})-g_{i}(x)\geq-
g_{i}(x),$

where we used the fact that $x^{*}$ is feasible for the last inequality. This
concludes that the summand in (23) is non-positive since
$\mathrm{d}\lambda_{i}\geq 0$. Thus, after some elementary manipulations, we
obtain the following upper bound on $\mathrm{d}W$:

$\displaystyle\mathrm{d}W\leq$ $\displaystyle-\alpha
W\mathrm{d}t-(2\delta-3\alpha/2)|u|^{2}\mathrm{d}t\leq-\alpha W\mathrm{d}t.$

Applying Grönwall’s inequality then implies the desired upper bound on
$f(x(t))-f(x^{*})$.

The lower bound is obtained from a perturbation analysis, using an argument
similar to Muehlebach and Jordan [26]. We define

$f^{*}(t):=\min_{\xi\in\mathbb{R}^{n}}f(\xi),\quad\text{s.t.}\quad
g(\xi)\geq\min\\{0,g(x(0))\\}e^{-\alpha t},$

which is of the form (1), but with a modified right-hand side of the
constraints. The trajectory $x(t)$ satisfies
$g(x(t))\geq\min\\{0,g(x(0))\\}e^{-\alpha t}$ and is therefore a feasible
candidate for the above minimization, which implies $f(x(t))\geq f^{*}(t)$.
The minimum is clearly attained, since $f$ is bounded below,
$f(x)\rightarrow\infty$ for $|x|\rightarrow\infty$, and the modified set of
feasible points is closed. The multiplier $\lambda^{*}$ satisfying the Karush-
Kuhn-Tucker conditions of (1) captures the sensitivity of the cost function
with respect to perturbations of the constraints. This means that
$-\lambda^{*}$ is guaranteed to satisfy the following inequality [see, e.g.,
40, p. 277]:

$f^{*}(t)-f^{*}\geq\min\\{0,g(x(0))\\}^{\mathsf{\tiny{T}}}\lambda^{*}e^{-\alpha
t},$

which combined with $f(x(t))\geq f^{*}(t)$ implies the desired lower bound.

#### ii) Nesterov - constant parameters ($\alpha=\sqrt{\mu}-\mu/2$):

We consider again the evolution of the function

$W(x,u)=\frac{1}{2}|a(x-x^{*})+u|^{2}+f(x),$

along the trajectories $(x(t),u(t))$, where $a:=\alpha$. The corresponding
differential measure is given by

$\mathrm{d}W=(a(x-x^{*})+(u^{+}+u^{-})/2)^{\mathsf{\tiny{T}}}\\\
(au\mathrm{d}t+\mathrm{d}u)+\nabla f(x)^{\mathsf{\tiny{T}}}u\mathrm{d}t.$

From the proof of Thm. 2 we obtain

$(u^{+}+u^{-})/2^{\mathsf{\tiny{T}}}\mathrm{d}u+\nabla
f(x)u\mathrm{d}t\leq-(2\delta+\mu\beta)|u|^{2}\mathrm{d}t\\\
+u^{\mathsf{\tiny{T}}}(\nabla f(x+\beta u)-\nabla f(x))\mathrm{d}t+\sum_{i\in
I_{x}}-\alpha g_{i}\mathrm{d}\lambda_{i}.$ (24)

Including the remaining terms yields the following upper bound on
$\mathrm{d}W$:

$-(2\delta-a)|u|^{2}\mathrm{d}t+(a^{2}-2\delta
a)(x-x^{*})^{\mathsf{\tiny{T}}}u\mathrm{d}t-(\nabla f(y)-\nabla
f(x))^{\mathsf{\tiny{T}}}u\mathrm{d}t-a(x-x^{*})^{\mathsf{\tiny{T}}}\nabla
f(y)\mathrm{d}t\\\ +\sum_{i\in I_{x}}(a\nabla
g_{i}(x)^{\mathsf{\tiny{T}}}(x-x^{*})-\alpha g_{i})\mathrm{d}\lambda_{i},$
(25)

where we introduced the variable $y:=x+\beta u$ to simplify notation. The fact
that $f$ is strongly convex means that the following holds:

$-(x-x^{*})^{\mathsf{\tiny{T}}}\nabla f(y)\leq-(f(x)-f(x^{*}))-\beta(\nabla
f(x)-\nabla f(y))^{\mathsf{\tiny{T}}}u-\mu\beta^{2}|u|^{2}\\\ -\mu\beta
u^{\mathsf{\tiny{T}}}(x-x^{*})-\frac{\mu}{2}|x-x^{*}|^{2}.$ (26)

In addition, $C$ is convex, which implies that each $g_{i}$ is concave. As a
result,

$\nabla g_{i}(x)^{\mathsf{\tiny{T}}}(x^{*}-x)\geq g_{i}(x^{*})-g_{i}(x)\geq-
g_{i}(x),$

where we used the fact that $x^{*}$ is feasible for the last inequality. The
summand in (25) can therefore be upper bounded by

$\left(a\nabla g_{i}(x)^{\mathsf{\tiny{T}}}(x-x^{*})-\alpha
g_{i}(x)\right)\mathrm{d}\lambda_{i}\leq(a-\alpha)g_{i}(x)\mathrm{d}\lambda_{i},$

since, by definition $\mathrm{d}\lambda_{i}\geq 0$. The upper bound vanishes
since $a=\alpha$. Combined with (26), and after some elementary manipulations,
this yields the following upper bound on $\mathrm{d}W$:

$\mathrm{d}W\leq-
aW\mathrm{d}t-(2\delta-3a/2+\mu\beta)|u|^{2}\mathrm{d}t-\frac{1}{2}(a\mu-a^{3})|x-x^{*}|^{2}\mathrm{d}t\\\
+(2a^{2}-2\delta a-a\mu\beta)(x-x^{*})^{\mathsf{\tiny{T}}}u\mathrm{d}t,$

where we have used the fact that $1-\beta a\geq 0$. We note that the term
$2a^{2}-2\delta a-a\mu\beta$ vanishes, that $a(\mu-a^{2})\geq 0$, and
$2\delta-3a/2+\mu\beta\geq 0$ for all $\mu\in[0,1]$. We therefore obtain
$\mathrm{d}W\leq-aW\mathrm{d}t$, which, by Grönwall’s inequality, implies the
desired result.

#### iii) Nesterov - varying parameters ($\alpha=2/(t+3)$):

The proof follows the same steps. We consider the evolution of the function

$\tilde{W}(x,u,t)=\frac{1}{2}|\tilde{a}(t)(x-x^{*})+u|^{2}+f(x),$

along the trajectories $(x(t),u(t))$, where we choose
$\tilde{a}(t):=\alpha(t)=2/(t+3)$. The differential measure
$\mathrm{d}\tilde{W}$ is given by

$\mathrm{d}\tilde{W}=(\tilde{a}(x-x^{*})+(u^{+}+u^{-})/2)^{\mathsf{\tiny{T}}}(\tilde{a}u\mathrm{d}t+\dot{\tilde{a}}(x-x^{*})\mathrm{d}t+\mathrm{d}u)+\nabla
f(x)^{\mathsf{\tiny{T}}}u\mathrm{d}t,$

where according to the chain rule the derivative of $\tilde{a}$ with respect
to $t$ enters. This leads to two additional quadratic terms of the type

$\tilde{a}\dot{\tilde{a}}|x-x^{*}|^{2}\mathrm{d}t\quad\text{and}\quad\dot{\tilde{a}}u^{\mathsf{\tiny{T}}}(x-x^{*})\mathrm{d}t.$

As a result, by following the same steps as in the variant ii) (see previous
section), we obtain

$\mathrm{d}\tilde{W}\leq-\tilde{a}\tilde{W}\mathrm{d}t-(2\delta-3\tilde{a}/2+\mu\beta)|u|^{2}\mathrm{d}t-\frac{1}{2}(\tilde{a}\mu-\tilde{a}^{3}-2\tilde{a}\dot{\tilde{a}})|x-x^{*}|^{2}\mathrm{d}t\\\
+(2\tilde{a}^{2}-2\delta\tilde{a}-\tilde{a}\mu\beta+\dot{\tilde{a}})(x-x^{*})^{\mathsf{\tiny{T}}}u\mathrm{d}t.$

We note that $2\tilde{a}^{2}-2\delta\tilde{a}+\dot{\tilde{a}}$ vanishes (in
fact $\alpha=\tilde{a}$ is deliberately constructed in this way). The same
applies to $2\delta-3\tilde{a}/2$ and
$\tilde{a}^{3}-2\tilde{a}\dot{\tilde{a}}$, which simplifies the above
inequality to

$\displaystyle\mathrm{d}\tilde{W}\leq-\tilde{a}\tilde{W}\mathrm{d}t-\mu\beta|u|^{2}\mathrm{d}t-\frac{1}{2}\tilde{a}\mu|x-x^{*}|^{2}\mathrm{d}t-\tilde{a}\mu\beta(x-x^{*})^{\mathsf{\tiny{T}}}u\mathrm{d}t.$

Applying Young’s inequality to the cross-term $(x-x^{*})^{\mathsf{\tiny{T}}}u$
concludes that $\mathrm{d}\tilde{W}\leq-\tilde{a}\tilde{W}\mathrm{d}t$ for any
$\mu\in[0,1]$. We finally apply Grönwall’s inequality, which yields

$\displaystyle\tilde{W}(x(t),u(t),t)\leq\tilde{W}(x(0),u(0),0)\exp\left(-\int_{0}^{t}\tilde{a}(s)\mathrm{d}s\right)=\tilde{W}(x(0),u(0),0)\frac{9}{(t+3)^{2}},$

and concludes the proof. ∎

##  Appendix F Convergence Analysis for (8)

We restrict ourselves to the case $\epsilon=0$, which allows us to restate the
algorithm in the following way:

$\displaystyle u_{k+1}$
$\displaystyle=\mathop{\mathrm{argmin}}_{v\in\mathbb{R}^{n}}\frac{1}{2}|v-u_{k}+2\delta
u_{k}T_{k}+\nabla f(x_{k}+\beta u_{k})T_{k}|^{2}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
\gamma_{i}(x_{k},v)\geq 0,\leavevmode\nobreak\ \leavevmode\nobreak\ i\in
I_{x_{k}},$ $\displaystyle x_{k+1}$ $\displaystyle=x_{k}+T_{k}u_{k+1}.$ (27)

We note that the following arguments still apply if the gradient is
stochastic, as long as the gradient perturbations have zero mean and bounded
second moments.

We start with the following three lemmas:

######  Lemma 8.

Let the assumptions of Thm. 4 be satisfied and let $\text{dist}(x,C)$ denote
the distance of $x\in\mathbb{R}^{n}$ to the set $C$, that is,
$\text{dist}(x,C)=\min_{y\in C}|y-x|$. Then, the iterates $x_{k}$ are bounded,
$\text{dist}(x_{k},C)\rightarrow 0$, and there exists a constant
$c_{\text{g}}>0$ and $k_{0}>0$ such that

$g_{i}(x_{k})\geq-c_{\text{g}}T_{k},$

for all $k\geq k_{0}$ and all $i\in\\{1,\dots,n_{\text{g}}\\}$.

###### Proof.

We start by considering the first constraint, $g_{1}$, since
$\min\\{g_{1}(x),0\\}$ has compact level sets. We distinguish the following
two cases.

i) $g_{1}(x_{k})<0$: As a consequence of smoothness we infer that
$g_{1}(x_{k+1})\geq g_{1}(x_{k})+T_{k}\nabla
g_{1}(x_{k})^{\mathsf{\tiny{T}}}u_{k+1}-T_{k}^{2}L_{\text{g}1}|u_{k+1}|^{2}/2$,
where $L_{\text{g}1}$ is the smoothness constant of $g_{1}$. Due to the fact
that $g_{1}(x_{k})<0$ we have we have $1\in I_{x_{k}}$ and $\nabla
g_{1}(x_{k})^{\mathsf{\tiny{T}}}u_{k+1}T_{k}\geq-\alpha T_{k}g_{1}(x_{k})$.
This means that

$g_{1}(x_{k+1})\geq g_{1}(x_{k})(1-\alpha
T_{k})-T_{k}^{2}L_{\text{g}1}c_{\text{u}}^{2}/2,$

where we used the upper bound $c_{\text{u}}$ on $u_{k+1}$. By the properties
of the sequence $T_{k}=T_{0}/k^{s}$ this implies that $g_{1}(x_{k+1})$ is
bounded below and

$g_{1}(x_{j})\geq-c_{\text{g}1}T_{j},$

for some constant $c_{\text{g}1}$ and $j=k+1,k+2,\dots$, as long as $1\in
I_{x_{j}}$, [see for example 41, Chapter 2, Lemma 4 and Lemma 5].

ii) $g_{1}(x_{k})\geq 0$ and $g_{1}(x_{k+1})<0$ (the constraint $g_{1}$
becomes active at $k+1$). We have again $g_{1}(x_{k+1})\geq\nabla
g_{1}(x_{k})^{\mathsf{\tiny{T}}}T_{k}u_{k+1}-T_{k}^{2}L_{\text{g}1}c_{\text{u}}^{2}/2$
as a result of the smoothness of $g_{1}$ and the boundedness of $u_{k+1}$. Due
to the fact that $g_{1}(x_{k})\geq 0$ we conclude that $x_{k}$ is bounded
($\min\\{g_{1}(x),0\\}$ has compact level sets), which means that $|\nabla
g_{1}(x_{k})|\leq\bar{c}_{\text{g1}}$ for some constant
$\bar{c}_{\text{g1}}>0$. Thus, we conclude

$g_{1}(x_{k+1})\geq-
T_{k}(\bar{c}_{\text{g1}}c_{\text{u}}+T_{k}L_{\text{g}1}c_{\text{u}}^{2}/2).$

By combining the two cases we infer that $g_{1}(x_{k})$ is bounded below for
all $k\geq 0$, which implies that $x_{k}$ is bounded. We can now apply a
similar reasoning to all the remaining constraints, which establishes that
$\text{dist}(x_{k},C)\rightarrow 0$ and $g_{i}(x_{k})\geq-c_{\text{g}}T_{k}$
for some constant $c_{\text{g}}>0$, all $i=1,2,\dots,n_{\text{g}}$, and all
$k$ large enough. ∎

######  Lemma 9.

Let $\bar{C}$ be a compact set. Then, there exists a constant $c_{\lambda}$
such that for all $x$ with $x\in\bar{C}$ and
$\lambda\in\mathbb{R}^{n_{\text{g}}}$ with $\lambda_{i}\geq 0$ $\forall i\in
I_{x}$ and $\lambda_{i}=0$ $\forall i\not\in I_{x}$,

$|\sum_{i=1}^{n_{\text{g}}}\nabla g_{i}(x)\lambda_{i}|\geq
c_{\lambda}|\lambda|.$

###### Proof.

The inequality is satisfied for $\lambda=0$. We therefore consider the case
$\lambda\neq 0$ and without loss of generality assume that $|\lambda|=1$. We
argue by contradiction and therefore assume that there are two convergent
sequences, $\lambda^{j}$ and $x_{j}$, with $|\lambda^{j}|=1$,
$\lambda^{j}_{i}\geq 0$ for all $i\in I_{x_{j}}$, $\lambda^{j}_{i}=0$ for all
$i\not\in I_{x_{j}}$, and $x_{j}\in C$, such that
$\sum_{i=1}^{n_{\text{g}}}\nabla g_{i}(x_{j})\lambda^{j}_{i}\rightarrow 0$.
Let the limit point of $x_{j}$ be denoted by $x$ and the limit point of
$\lambda^{j}$ by $\lambda$. As a consequence of constraint qualification,
there exists a vector $w\in\mathbb{R}^{n}$ such that
$w^{\mathsf{\tiny{T}}}\nabla g_{i}(x)\lambda_{i}>0$ for all $i\in I_{x}$. We
note that $g_{i}(x)>0$ implies $g_{i}(x_{j})>0$ for large enough $j$. This
shows that $I_{x_{j}}\subset I_{x}$ and therefore

$\sum_{i\in I_{x_{j}}}w^{\mathsf{\tiny{T}}}\nabla g_{i}(x)\lambda_{i}>c_{1}$

for all large enough $j$ and a small enough $c_{1}>0$. However, by continuity
of $\nabla g$ we have

$\sum_{i=1}^{n_{\text{g}}}w^{\mathsf{\tiny{T}}}\nabla
g_{i}(x_{j})\lambda_{i}^{j}=\sum_{i\in I_{x_{j}}}w^{\mathsf{\tiny{T}}}\nabla
g_{i}(x)\lambda_{i}+\sum_{i\in I_{x_{j}}}w^{\mathsf{\tiny{T}}}\nabla
g_{i}(x)(\lambda_{i}^{j}-\lambda_{i})\\\ +\sum_{i\in
I_{x_{j}}}w^{\mathsf{\tiny{T}}}(\nabla g_{i}(x_{j})-\nabla
g_{i}(x))\lambda_{i}^{j}>c_{1}/2,$

for all large enough $j$, which is a contradiction. ∎

######  Lemma 10.

Let the assumptions of Thm. 4 be satisfied. If $x_{k(j)}\rightarrow\bar{x}\in
C$, $u_{k(j)}\rightarrow 0$, and $R_{k(j)}/T_{k(j)}\rightarrow\bar{R}$ for a
subsequence $k(j),j=1,2,\dots$, then $\bar{x}$ and $\bar{R}$ satisfy

$-\nabla f(\bar{x})+\bar{R}=0,\quad-\bar{R}\in N_{V_{\alpha}(\bar{x})}(0).$

###### Proof.

We conclude from the continuity of $g_{i}$ that $g_{i}(\bar{x})>0$ implies
that $g_{i}(x_{k(j)})>0$ for all $j$ large enough and all
$i\in\\{1,2,\dots,n_{\text{g}}\\}$. This means that $I_{x_{k(j)}}\subset
I_{\bar{x}}$ for all $j$ large enough. We define a slightly modified version
of $V_{\alpha}(x)$ as follows

$\tilde{V}_{\alpha}(x):=\\{u\in\mathbb{R}^{n}\leavevmode\nobreak\
|\leavevmode\nobreak\ \nabla g_{i}(x)^{\mathsf{\tiny{T}}}u\geq-\alpha
g_{i}(x),\quad\forall i\in I_{\bar{x}}\\},$

which ensures that $V_{\alpha}(x_{k(j)})\supset\tilde{V}_{\alpha}(x_{k(j)})$
for large $j$. Hence, the corresponding normal cones satisfy

$N_{V_{\alpha}(x_{k(j)})}(u)\subset N_{\tilde{V}_{\alpha}(x_{k(j)})}(u),$

for all $u\in\tilde{V}_{\alpha}(x_{k(j)})$, which implies, by the update rule
of algorithm (8),

$-\frac{R_{k(j)}}{T_{k(j)}}\in N_{V_{\alpha}(x_{k(j)})}(u_{k(j)+1})\subset
N_{\tilde{V}_{\alpha}(x_{k(j)})}(u_{k(j)+1}),$

where we have used the fact that the normal cone is a cone. We now show that
this implies $-\bar{R}\in N_{\tilde{V}_{\alpha}(\bar{x})}(0)$ and argue by
contradiction. This means that there exists a
$\hat{u}\in\tilde{V}_{\alpha}(\bar{x})$ such that
$-\bar{R}^{\mathsf{\tiny{T}}}\hat{u}>c_{1}$ for a small $c_{1}>0$. From
constraint qualification, we infer that there exists a $w\in\mathbb{R}^{n}$
and an $\varepsilon>0$ such that $\hat{u}+\varepsilon
w\in\tilde{V}_{\alpha}(x_{k(j)})$ for all $j$ sufficiently large with
$-\bar{R}^{\mathsf{\tiny{T}}}(\hat{u}+\varepsilon w)>c_{1}/2$. However, this
leads to a contradiction, since

$0\geq-\frac{R_{k(j)}}{T_{k(j)}}^{\mathsf{\tiny{T}}}(\hat{u}+\varepsilon
w-u_{k(j)+1})=\\\ \underbrace{-\bar{R}^{\mathsf{\tiny{T}}}(\hat{u}+\varepsilon
w)}_{>c_{1}/2>0}+\underbrace{\bar{R}^{\mathsf{\tiny{T}}}u_{k(j)+1}-\left(\frac{R_{k(j)}}{T_{k(j)}}-\bar{R}\right)^{\mathsf{\tiny{T}}}(\hat{u}+\varepsilon
w-u_{k(j)+1})}_{\rightarrow 0}.$

This shows that $-\bar{R}\in N_{\tilde{V}_{\alpha}(\bar{x})}(0)$ and the
desired result follows from the fact that
$\tilde{V}_{\alpha}(\bar{x})=V_{\alpha}(\bar{x})$. ∎

The following proposition will be based on dividing the update (8) into the
following two steps

$\displaystyle\bar{x}_{k}$
$\displaystyle=x_{k},\qquad\qquad\qquad\qquad\qquad\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
x_{k+1}=\bar{x}_{k}+u_{k+1}T_{k},$ $\displaystyle\bar{u}_{k}$
$\displaystyle=u_{k}-T_{k}f_{\text{d}}(x_{k},u_{k})+R_{k},\qquad
u_{k+1}=\bar{u}_{k}-T_{k}\nabla f(\bar{x}_{k}),$ (28)

where $f_{\text{d}}(x_{k},u_{k}):=2\delta u_{k}+\nabla f(x_{k}+\beta
u_{k})-\nabla f(x_{k})$ contains the dissipative terms. The first step, which
maps $(x_{k},u_{k})$ to $(\bar{x}_{k},\bar{u}_{k})$ is an update of the
velocity with the dissipative terms and the constraint forces $R_{k}$, whereas
the second step, which maps $(\bar{x}_{k},\bar{u}_{k})$ to $(x_{k+1},u_{k+1})$
is a symplectic Euler discretization that captures the conservative parts of
the underlying dynamics. In the following we will exploit the fact that the
second step is a symplectic map. We refer the reader to [9] and [7] for
additional details.

We now prove the following lemma, which will lie at the heart of the
convergence proof in discrete time.

######  Lemma 11.

Let the assumptions of Thm. 4 be satisfied. Then, there exists a function
$V_{k}(x,u)$, which is bounded below (uniformly in $k$), such that

$V_{k+1}(x_{k+1},u_{k+1})-V_{k}(x_{k},u_{k})\leq-
c_{\text{V}1}T_{k}|u_{k}|^{2}-c_{\text{V}2}|R_{k}-T_{k}\nabla
f(x_{k})|^{2}-c_{\text{V}3}T_{k}|R_{k}|^{2}\\\
+c_{\text{V}4}T_{k}^{2}|u_{k+1}|+c_{\text{V}5}T_{k}^{3}-\alpha\sum_{i\in
I_{x_{k}}}\lambda_{k}^{i}g_{i}(x_{k}),$

for all $k$ large enough, where
$c_{\text{V}1},c_{\text{V}2},c_{\text{V}3},c_{\text{V}4},c_{\text{V}5}>0$ are
constant.

###### Proof.

From Lemma 8 we infer that $x_{k}$ is bounded and therefore contained in a
compact set, which we denote by $\bar{C}$. Without loss of generality we
assume that $f$ is analytic on $\bar{C}$ (note that by the Stone-Weierstrass
theorem we can find a polynomial that approximates $f$ and $\nabla f$
arbitrarily closely, see [9] for details). We can now invoke Muehlebach and
Jordan [Proposition 9 of 9], which constructs the function $F_{k}(x,u)$
(continuous in $x$ and $u$) such that

$|\tilde{H}_{k}(x_{k+1},u_{k+1})-\tilde{H}_{k}(\bar{x}_{k},\bar{u}_{k})|\leq
c_{\text{E}}T_{k}^{4}$ (29)

for $k$ sufficiently large, where

$\tilde{H}_{k}(x,u):=\frac{1}{2}|u|^{2}+f(x)-\frac{T_{k}}{2}\nabla
f(x)^{\mathsf{\tiny{T}}}u+T_{k}^{2}F_{k}(x,u),$

and $c_{\text{E}}$ is constant. The function $F_{k}(x,u)$ is guaranteed to be
bounded for $x\in\bar{C}$, $|u|\leq c_{\text{u}}$ (uniformly in $k$). We note
that (29) implies that $\tilde{H}$ is almost conserved when applying the
symplectic step in (28) (in fact, the right hand side can be replaced with
$(T_{k})^{q}$ for an arbitrarily large $q>0$). The function $\tilde{H}$ is
therefore often referred to as a modified Hamiltonian and arises from the fact
that the second step in (28) is a symplectic map.

We claim that the function

$V_{k}(x,u):=\tilde{H}_{k}(x,u)+\frac{T_{k}d_{k}}{2}\nabla
f(x)^{\mathsf{\tiny{T}}}u$

satisfies the desired properties, where $d_{k}:=1-\bar{d}T_{k}>0$ ($\bar{d}>0$
is bounded and will be chosen subsequently). The choice $d_{k}=1-\bar{d}T_{k}$
is motivated by the fact that, as a result, $V_{k+1}(x,u)-V_{k}(x,u)$ is of
order $T_{k}^{2}$ for $x_{k}\in\bar{C}$ and $|u_{k}|\leq c_{\text{u}}$. We
start by considering step 1 in (28), where we note that

$\displaystyle\frac{1}{2}|\bar{u}_{k}|^{2}$
$\displaystyle-\frac{1}{2}|u_{k}|^{2}=\frac{1}{2}(\bar{u}_{k}-u_{k})^{\mathsf{\tiny{T}}}(\bar{u}_{k}+u_{k})$
$\displaystyle=-T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}u_{k}+\frac{1}{2}u_{k}^{\mathsf{\tiny{T}}}R_{k}+\frac{1}{2}T_{k}^{2}|f_{\text{d}}|^{2}-\frac{3}{4}T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}R_{k}+\frac{1}{4}|R_{k}|^{2}$
$\displaystyle=-T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}u_{k}+\frac{1}{2}u_{k+1}^{\mathsf{\tiny{T}}}R_{k}+\frac{1}{2}T_{k}^{2}|f_{\text{d}}|^{2}-\frac{1}{4}T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}R_{k}-\frac{1}{4}|R_{k}|^{2}+\frac{1}{2}T_{k}R_{k}^{\mathsf{\tiny{T}}}\nabla
f(x_{k}),$

where we replaced $u_{k}$ by $u_{k+1}$ in the last step and omitted the
arguments of $f_{\text{d}}(x_{k},u_{k})$. We further have

$\displaystyle-\frac{T_{k}}{2}(1-d_{k})\nabla
f(x_{k})^{\mathsf{\tiny{T}}}(\bar{u}_{k}-u_{k})$
$\displaystyle=-\frac{T_{k}}{2}(1-d_{k})\nabla
f(x_{k})^{\mathsf{\tiny{T}}}(-T_{k}f_{\text{d}}+R_{k}),$

which means that over step 1 in (28) the function $V_{k}$ changes by

$\displaystyle V_{k}(\bar{x}_{k},\bar{u}_{k})-V_{k}(x_{k},u_{k})$
$\displaystyle\leq-
T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}u_{k}+\frac{1}{2}u_{k+1}^{\mathsf{\tiny{T}}}R_{k}+\frac{1}{2}T_{k}^{2}|f_{\text{d}}|^{2}-\frac{1}{4}T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}R_{k}-\frac{1}{4}|R_{k}|^{2}$
$\displaystyle+\frac{d_{k}}{2}T_{k}R_{k}^{\mathsf{\tiny{T}}}\nabla
f(x_{k})+\frac{T_{k}^{2}(1-d_{k})}{2}\nabla
f(x_{k})^{\mathsf{\tiny{T}}}f_{\text{d}}$
$\displaystyle+T_{k}^{2}\nabla_{u}F_{k}(x_{k},u_{k})^{\mathsf{\tiny{T}}}(-T_{k}f_{\text{d}}+R_{k})+T_{k}^{2}c_{\text{F}}|-T_{k}f_{\text{d}}+R_{k}|^{2},$

where $c_{\text{F}}>0$ is a bound on the smoothness constant of $F_{k}$. We
now proceed to the second step of (28), where we exploit the fact that
$\tilde{H}_{k}$ is invariant up to terms of order $T_{k}^{4}$. We are
therefore left with analyzing the term $T_{k}d_{k}\nabla
f(x)^{\mathsf{\tiny{T}}}u$, which gives

$\displaystyle V_{k}(x_{k+1},u_{k+1})\\!-\\!V_{k}(\bar{x}_{k},\bar{u}_{k})$
$\displaystyle\\!\leq\\!\frac{T_{k}d_{k}}{2}(\nabla f(x_{k+1})\\!-\\!\nabla
f(x_{k}))^{\mathsf{\tiny{T}}}\\!u_{k+1}\\!-\\!\frac{T_{k}^{2}d_{1}}{2}|\nabla
f(x_{k})|^{2}\\!+\\!c_{\text{E}}T_{k}^{4}$
$\displaystyle\leq\frac{T_{k}^{2}d_{k}}{2}|u_{k+1}|^{2}-\frac{T_{k}^{2}d_{k}}{2}|\nabla
f(x_{k})|^{2}+c_{\text{E}}T_{k}^{4},$

where we used the fact that $\nabla f$ is 1-smooth in the second step. We
therefore obtain

$V_{k}(x_{k+1},u_{k+1})-V_{k}(x_{k},u_{k})\leq-
T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}u_{k}+\frac{1}{2}u_{k+1}^{\mathsf{\tiny{T}}}R_{k}-\frac{1}{4}T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}R_{k}+\frac{1}{2}T_{k}^{2}|f_{\text{d}}|^{2}\\\
+\frac{T_{k}^{2}d_{k}}{2}|u_{k+1}|^{2}-\frac{T_{k}\bar{d}}{4}|R_{k}|^{2}-\frac{d_{k}}{4}|R_{k}-T_{k}\nabla
f(x_{k})|^{2}-\frac{d_{k}}{4}T_{k}^{2}|\nabla f(x_{k})|^{2}\\\
+c_{\text{E}4}T_{k}^{4}+T_{k}^{2}\nabla_{u}F_{k}(x_{k},u_{k})^{\mathsf{\tiny{T}}}(-T_{k}f_{\text{d}}+R_{k})+T_{k}^{2}c_{\text{F}}|-T_{k}f_{\text{d}}+R_{k}|^{2}.$

In addition we can bound $f_{\text{d}}^{\mathsf{\tiny{T}}}u_{k}$ and
$|f_{\text{d}}|^{2}$ as follows:

$\displaystyle-
f_{\text{d}}^{\mathsf{\tiny{T}}}u_{k}\leq-\bar{\delta}|u_{k}|^{2},\quad|f_{\text{d}}|^{2}\leq(2\delta+\beta)^{2}|u_{k}|^{2},$

where $\bar{\delta}=2\delta$ in the convex case and
$\bar{\delta}=2\delta-\beta$ in the nonconvex case. We can thus apply Young’s
inequality to the term $f_{\text{d}}^{\mathsf{\tiny{T}}}R_{k}$, which yields

$\displaystyle-\frac{1}{4}T_{k}f_{\text{d}}^{\mathsf{\tiny{T}}}R_{k}\leq\frac{1}{2}T_{k}\bar{\delta}|u_{k}|^{2}+T_{k}\frac{|2\delta+\beta|^{2}}{32\bar{\delta}}|R_{k}|^{2}.$

We further note that $|\nabla_{u}F_{k}(x_{k},u_{k})|\leq
c_{\text{F}2}(|u_{k}|+|\nabla f(x_{k})|)$ and apply Young’s inequality to all
cross-terms. This leads to

$\displaystyle
V_{k}(x_{k+1},u_{k+1})-V_{k}(x_{k},u_{k})\leq-\frac{1}{2}T_{k}\bar{\delta}|u_{k}|^{2}-\frac{T_{k}\bar{d}}{8}|R_{k}|^{2}-\frac{d_{k}}{4}|R_{k}-T_{k}\nabla
f(x_{k})|^{2}$
$\displaystyle+\frac{T_{k}^{2}d_{k}}{2}|u_{k+1}|^{2}+\frac{1}{2}u_{k+1}^{\mathsf{\tiny{T}}}R_{k}+c_{\text{E}}T_{k}^{4},$

for large enough $k$ and for a constant $\bar{d}$ large enough. In addition,
the change of $V_{k+1}(x_{k+1},u_{k+1})-V_{k}(x_{k+1},u_{k+1})$ can be shown
to be of the order $\mathcal{O}(T_{k}^{2}|u_{k+1}|+T_{k}^{3})$. We further
infer from the update equation for $u_{k+1}$ that
$u_{k+1}^{\mathsf{\tiny{T}}}R_{k}=-\alpha\sum_{i\in
I_{x_{k}}}\lambda_{k}^{i}g_{i}(x_{k})$, which yields the desired result. ∎

######  Proposition 12.

Let the assumptions of Thm. 4 be satisfied. Then $x_{k}$ converges to the set
of stationary points of (1) and $u_{k}$ converges to zero.

###### Proof.

From Lemma 11 we infer the existence of a function $V_{k}(x,u)$ (bounded
below, uniformly in $k$) such that

$V_{k+1}(x_{k+1},u_{k+1})-V_{k}(x_{k},u_{k})\leq-
c_{\text{V}1}T_{k}|u_{k}|^{2}-c_{\text{V}2}|R_{k}|^{2}+c_{\text{V}3}T_{k}^{2}-\alpha\sum_{i\in
I_{x_{k}}}\lambda_{k}^{i}g_{i}(x_{k})$

for $k\geq k_{0}$ and for some constants $k_{0}>0$,
$c_{\text{V}1},c_{\text{V}2},c_{\text{V}3}>0$, where we have used the fact
that $x_{k}$ and $u_{k}$ are bounded. (We have abused notation by reusing
$c_{\text{V}1},c_{\text{V}2},c_{\text{V}3}$.) We now invoke Lemma 8 and Lemma
9 for bounding the terms $\lambda_{k}^{i}g_{i}(x_{k})$

$|-\alpha\sum_{i\in
I_{x_{k}}}\lambda_{k}^{i}g_{i}(x_{k})|\leq\frac{\alpha}{c_{\lambda}}c_{\text{g}}T_{k}|R_{k}|n_{\text{g}}\leq\frac{c_{\text{V}2}}{2}|R_{k}|^{2}+\frac{1}{2c_{\text{V}2}}(n_{\text{g}}c_{\text{g}}\frac{\alpha}{c_{\lambda}})^{2}T_{k}^{2},$
(30)

where we have used Young’s inequality in the second step. We therefore obtain

$V_{k+1}(x_{k+1},u_{k+1})-V_{k}(x_{k},u_{k})\leq-
c_{\text{V}1}T_{k}|u_{k}|^{2}-\frac{c_{\text{V}2}}{2}|R_{k}|^{2}+\bar{c}_{\text{V}3}T_{k}^{2},$
(31)

for large enough $k$ and a modified constant $\bar{c}_{\text{V}3}>0$. We
further introduce the constants

$-c_{\text{V}0}=\inf_{k}V_{k}(x_{k},u_{k}),\qquad
c_{\text{T}1}=\sum_{k=k_{0}}^{\infty}T_{k}^{2},$

and make the following claim: There exists a subsequence $k(j),j=1,2,\dots$,
such that

$-c_{\text{V}1}T_{k(j)}|u_{k(j)}|^{2}-\frac{c_{\text{V}2}}{2}|R_{k(j)}|^{2}\geq\underbrace{\left(-\bar{c}_{\text{V}3}-\frac{V_{k_{0}}(x_{k_{0}},u_{k_{0}})+c_{\text{V}0}}{c_{\text{T}1}}-\frac{1}{c_{\text{T}1}}\right)}_{:=-c_{\text{V}4}}T_{k(j)}^{2},$

where $c_{\text{V}4}>0$. For the sake of contradiction we assume that the
claim is not true, which means that

$-c_{\text{V}1}T_{k}|u_{k}|^{2}-\frac{c_{\text{V}2}}{2}|R_{k}|^{2}\leq-{c}_{\text{V}4}T_{k}^{2}$

for all $k>0$. However, we have chosen the constant $c_{\text{V}4}$
deliberately in such a way that we can generate a contradiction when summing
over $V_{k+1}(x_{k+1},u_{k+1})-V_{k}(x_{k},u_{k})$. More precisely,

$\displaystyle-c_{\text{V}0}-V_{k_{0}}(x_{k_{0}},u_{k_{0}})$
$\displaystyle\leq
V_{N}(x_{N},u_{N})-V_{k_{0}}(x_{k_{0}},u_{k_{0}})\leq\sum_{k=k_{0}}^{N}(-c_{\text{V}4}+\bar{c}_{\text{V}3})T_{k}^{2}$
$\displaystyle\leq\frac{-c_{\text{V}0}-V_{k_{0}}(x_{k_{0}},u_{k_{0}})-1}{c_{\text{T}1}}\sum_{k=k_{0}}^{N}T_{k}^{2}.$

which leads to the desired contradiction since the right-hand side approaches
$-c_{\text{V}_{0}}-V_{k_{0}}(x_{k_{0}},u_{k_{0}})-1$ for $N\rightarrow\infty$.

Upon passing to another subsequence we infer that $u_{k(j)}\rightarrow 0$,
$R_{k(j)}/T_{k(j)}\rightarrow\bar{R}$, and $x_{k(j)}\rightarrow\bar{x}\in C$.
By invoking Lemma 10 we conclude that $\bar{x}$ is stationary and $\bar{x}$
and $\bar{R}$ satisfy the Karush-Kuhn-Tucker conditions.

We now prove that the entire sequence converges. We infer from (31) that for
any $\varepsilon>0$ we can find an $m>0$ such that $V_{N}(x_{N},u_{N})\leq
V_{m}(x_{m},u_{m})+\varepsilon$ for all $N>m$. More precisely,

$\sum_{k=m}^{N-1}V_{k+1}(x_{k+1},u_{k+1})-V_{k}(x_{k},u_{k})\leq\bar{c}_{\text{V}3}T_{m}^{q}\sum_{k=m}^{N-1}\left(\frac{T_{k}}{T_{m}}\right)^{q}T_{k}^{2-q}\leq\bar{c}_{\text{V}3}T_{m}^{q}\sum_{k=m}^{\infty}T_{k}^{2-q},$
(32)

where the right-hand side is of the order $T_{m}^{q}$, since $T_{k}^{2-q}$ is
summable for a small enough $q>0$. This implies that all accumulation points
of $x_{k}$, $u_{k}$ must be contained in the same level set of the function
$|u|^{2}/2+f(x)=\lim_{k\rightarrow\infty}V_{k}(x,u)$. Thus, the sequence
$V_{k}(x_{k},u_{k})$ converges and satisfies

$\lim_{k\rightarrow\infty}V_{k}(x_{k},u_{k})=\lim_{k\rightarrow\infty}\left(\frac{1}{2}|u_{k}|^{2}+f(x_{k})\right)=f(\bar{x}),$

where $\bar{x}$ is the limit of $x_{k(j)}$ (see above). We further introduce
the collection of all accumulation points of $x_{k},u_{k}$, that is,

$\omega:=\bigcap_{m\geq 0}\text{cl}\\{(x_{k},u_{k}):k>m\\},$

where cl denotes closure. We note that $\omega\subset C\times\mathbb{R}^{n}$
is connected (since $|x_{k+1}-x_{k}|\rightarrow 0$ and
$|u_{k+1}-u_{k}|\rightarrow 0$) and for any $(x,u)\in\omega$,
$|u|^{2}/2+f(x)=f(\bar{x})$.

We claim that $\omega=\\{(\bar{x},0)\\}$. For the sake of contradiction, we
assume the existence of a sequence
$(\tilde{x}_{k},\tilde{u}_{k})\rightarrow(\bar{x},0)$,
$(\tilde{x}_{k},\tilde{u}_{k})\in\omega$ with
$(\tilde{x}_{k},\tilde{u}_{k})\neq(\bar{x},0)$ ($\omega$ is connected). We
consider first the case where $\bar{x}$ lies in the interior of $C$
(constraints are not active). In that case $\tilde{x}_{k}$ lies likewise in
the interior of $C$ for large $k$. However, in the absence of constraints the
dynamics in (8) are smooth, which means that $\tilde{x}_{k}$ and
$\tilde{u}_{k}$ are guaranteed to be equilibria and therefore satisfy
$\tilde{u}_{k}=0$ and $\nabla f(\tilde{x}_{k})=0$. However, this contradicts
the fact that $f$ has isolated stationary points. Next, we consider the case
where $\bar{x}$ lies on the boundary of $C$, where we infer from
$(\tilde{x}_{k},\tilde{u}_{k})\in\omega$ that
$f(\tilde{x}_{k})=f(\bar{x})-|\tilde{u}_{k}|^{2}/2$. Therefore
$f(\tilde{x}_{k})\leq f(\bar{x})$, $\tilde{x}_{k}\in C$ for all $k>0$, which
also contradicts the fact that $\bar{x}$ is an isolated stationary point.

Thus, we conclude $\omega=\\{(\bar{x},0)\\}$ and the result follows. ∎

##  Appendix G Proof of Thm. 5

We simplify the presentation by assuming that $L_{l}=T=1$, which is without
loss of generality. The proof hinges on the following Lyapunov function:

$V(x_{k},u_{k}):=\frac{1}{2}|\delta(x_{k}-x^{*})+(1-\delta)u_{k}|^{2}+l(x_{k})-l(x^{*}),$

which can be readily verified to be positive definite for all $\kappa_{l}\geq
1$. Moreover, the evolution of $V$ along the iterates (12) is given by

$\displaystyle V_{k+1}-V_{k}\leq-\delta(1-\frac{3\delta}{2})|u_{k}|^{2}-\beta
u_{k}^{\mathsf{\tiny{T}}}\nabla\bar{l}(y_{k})+\delta(x^{*}-x_{k})^{\mathsf{\tiny{T}}}\nabla\bar{l}(y_{k})$
$\displaystyle+\delta^{2}u_{k}^{\mathsf{\tiny{T}}}(x^{*}-x_{k})+\frac{1}{2}|\nabla\bar{l}(y_{k})|^{2}+l(x_{k+1})-l(x_{k}),$
(33)

where we abbreviated $V(x_{k},u_{k})$ by $V_{k}$. We have further used the
fact that $x_{k+1}-y_{k}=-\nabla f(y_{k})+\nabla
g(y_{k})\lambda_{k}:=-\nabla\bar{l}(y_{k})$, with $\lambda_{k}$ a set of
optimal multipliers of (12) (see also (34) below).

Next, we will relate $l(x_{k+1})$ to $l(y_{k})$. To do so, we start by
slightly reformulating the minimization in (12) by performing a change of
variables, $x=x_{k}+Tv$, which yields

$\displaystyle x_{k+1}$
$\displaystyle=\mathop{\mathrm{argmin}}_{x\in\mathbb{R}^{n}}\frac{1}{2}|x-y_{k}+\nabla
f(y_{k})/L_{l}|^{2}$ (34) $\displaystyle\text{s.t.}\leavevmode\nobreak\ \nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}(x-y_{k})\geq(1-\alpha
T)(g_{i}(x_{k})-g_{i}(y_{k}))-\alpha Tg_{i}(y_{k}),\leavevmode\nobreak\
\leavevmode\nobreak\ i\in[n_{\text{g}}],$

where we have used the fact that $T(1-2\delta T)=\beta$ to simplify the
objective function. We note that $x=\alpha(x^{*}-x_{k})+x_{k}$ is a feasible
solution candidate in the minimization (34). This follows from the following
reasoning:

$\displaystyle\nabla g_{i}(y_{k})(\alpha(x^{*}-x_{k})+x_{k}-y_{k})$
$\displaystyle=(1-\alpha)\nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}(x_{k}-y_{k})+\alpha\nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}(x^{*}-y_{k})$
$\displaystyle\geq(1-\alpha)(g_{i}(x_{k})-g_{i}(y_{k}))+\alpha(g_{i}(x^{*})-g(y_{k})),$

where concavity of $g$ is used in the second step together with the fact that
$T=1$ and $g_{i}(x^{*})\geq 0$. As a result, the stationarity condition of
(34) implies

$(x_{k+1}-y_{k}+\nabla
f(y_{k}))^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})+x_{k}-x_{k+1})\geq 0.$

This inequality can be restated as follows:

$(x_{k+1}-y_{k}+\nabla f(y_{k}))^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k}-(x_{k+1}-y_{k}))\\\
=-|\nabla\bar{l}(y_{k})|^{2}-\nabla\bar{l}(y_{k})^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k})\\\ +\nabla f(y_{k})^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k})-\nabla f(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k})\geq 0.$

Rearranging this inequality yields the following bound:

$\nabla
l(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k})\leq-|\nabla\bar{l}(y_{k})|^{2}-\nabla\bar{l}(y_{k})^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k})\\\ +\nabla f(y_{k})^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k})-{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k}).$ (35)

As a result of this inequality and the smoothness of $l$, we can bound the
evolution of $l(x_{k+1})$ as

$\displaystyle l(x_{k+1})$ $\displaystyle\leq l(y_{k})+\nabla
l(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k})+\frac{1}{2}|\nabla\bar{l}(y_{k})|^{2}$
$\displaystyle\leq
l(y_{k})-\frac{1}{2}|\nabla\bar{l}(y_{k})|^{2}-\nabla\bar{l}(y_{k})^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k})$ $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ +\nabla
f(y_{k})^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k})-{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k}).$ (36)

By combining (36) with (33) we obtain

$\displaystyle V_{k+1}-V_{k}\leq-\delta(1-\frac{3\delta}{2})|u_{k}|^{2}+\nabla
f(y_{k})^{\mathsf{\tiny{T}}}(\alpha(x^{*}-x_{k})-\beta
u_{k})+l(y_{k})-l(x_{k})$
$\displaystyle-{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k})+\delta^{2}u_{k}^{\mathsf{\tiny{T}}}(x^{*}-x_{k}),$
(37)

where we have exploited that $\alpha=\delta$ and the fact that all the terms
containing $\nabla\bar{l}(y_{k})$ cancel out (in fact the Lyapunov function
and the algorithm (12) are specifically engineered in this way). The strong
convexity of $f$ implies that

$\nabla f(y_{k})^{\mathsf{\tiny{T}}}(\delta(x^{*}-x_{k})-\beta
u_{k})\leq-\delta(f(x_{k})-f(x^{*}))-\frac{\delta\mu}{2}|x^{*}-x_{k}|^{2}-\delta\beta\mu(x_{k}-x^{*})^{\mathsf{\tiny{T}}}u_{k}\\\
-\frac{\mu}{2}\beta^{2}|u_{k}|^{2}+f(x_{k})-f(y_{k});$ (38)

see Muehlebach and Jordan [7, App. A6], which yields, combined with (37),

$V_{k+1}-V_{k}\leq-\delta(1-\frac{3\delta}{2}+\frac{\beta^{2}\mu}{2\delta})|u_{k}|^{2}-\frac{\delta\mu}{2}|x^{*}-x_{k}|^{2}+\delta(\beta\mu+\delta)(x^{*}-x_{k})^{\mathsf{\tiny{T}}}u_{k}\\\
-\delta(f(x_{k})-f(x^{*}))-{\lambda^{*}}g(y_{k})+{\lambda^{*}}^{\mathsf{\tiny{T}}}g(x_{k})-{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k}).$

This can be rearranged as follows:

$\\!\\!\\!V_{k+1}-V_{k}\leq-\delta
V_{k}-\frac{1}{2}(\delta-\delta^{2}-\delta^{3}+\beta^{2}\mu)|u_{k}|^{2}-\frac{\delta\mu-\delta^{3}}{2}|x^{*}-x_{k}|^{2}+\delta(\beta\mu+\delta^{2})(x^{*}-x_{k})^{\mathsf{\tiny{T}}}u_{k}\\\
+(1-\delta){\lambda^{*}}^{\mathsf{\tiny{T}}}g(x_{k})-{\lambda^{*}}^{\mathsf{\tiny{T}}}g(y_{k})-{\lambda^{*}}^{\mathsf{\tiny{T}}}\nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k}).$

We further note that the inequality constraint in (34) implies

$-\nabla
g_{i}(y_{k})^{\mathsf{\tiny{T}}}(x_{k+1}-y_{k})+(1-\delta)g(x_{k})-g(y_{k})\leq
0,$

and therefore we obtain

$V_{k+1}-V_{k}\leq-\delta
V_{k}-\frac{1}{2}(\delta-\delta^{2}-\delta^{3}+\beta^{2}\mu)|u_{k}|^{2}-\frac{\delta\mu-\delta^{3}}{2}|x^{*}-x_{k}|^{2}+\delta(\beta\mu+\delta^{2})(x^{*}-x_{k})^{\mathsf{\tiny{T}}}u_{k}.$

We note that $\beta\mu=(1-\mu)\delta^{2}$ and therefore
$\beta\mu+\delta^{2}=\delta^{2}(2-\mu)$. We now apply Young’s inequality to
the term $(x^{*}-x_{k})^{\mathsf{\tiny{T}}}u_{k}$ and conclude

$-\frac{1}{2}(\delta-\delta^{2}-\delta^{3}+\beta^{2}\mu)|u_{k}|^{2}-\frac{\delta\mu-\delta^{3}}{2}|x^{*}-x_{k}|^{2}+\delta^{3}(2-\mu)(x^{*}-x_{k})^{\mathsf{\tiny{T}}}u_{k}\leq
0,$

for all $\kappa_{l}\geq 1$, $x^{*}-x_{k}$, and $u_{k}$. This yields
$V_{k+1}-V_{k}\leq-\delta V_{k}$ and concludes the proof. ∎

## References

  *   * Su et al. [2016] Su, W., Boyd, S., Candès, E.J.: A differential equation for modeling Nesterov’s accelerated gradient method: Theory and insights. Journal of Machine Learning Research 17(153), 1–43 (2016)
  * Wibisono et al. [2016] Wibisono, A., Wilson, A.C., Jordan, M.I.: A variational perspective on accelerated methods in optimization. Proceedings of the National Academy of Sciences 113(47), 7351–7358 (2016)
  * Diakonikolas and Jordan [2021] Diakonikolas, J., Jordan, M.I.: Generalized momentum-based methods: A Hamiltonian perspective. SIAM Journal on Optimization 31(1), 915–944 (2021)
  * Krichene et al. [2015] Krichene, W., Bayen, A.M., Bartlett, P.L.: Accelerated mirror descent in continuous and discrete time. Advances in Neural Information Processing Systems 28, 2845–2853 (2015)
  * França et al. [2020] França, G., Sulam, J., Robinson, D.P., Vidal, R.: Conformal symplectic and relativistic optimization. Journal of Statistical Mechanics: Theory and Experiment 2020(12), 1–30 (2020)
  * Betancourt et al. [2018] Betancourt, M., Jordan, M.I., Wilson, A.C.: On symplectic optimization. arXiv:1802.03653v2, 1–20 (2018)
  * Muehlebach and Jordan [2019] Muehlebach, M., Jordan, M.I.: A dynamical systems perspective on Nesterov acceleration. Proceedings of Machine Learning Research 97, 4656–4662 (2019)
  * Muehlebach and Jordan [2020] Muehlebach, M., Jordan, M.I.: Continuous-time lower bounds for gradient-based algorithms. Proceedings of Machine Learning Research 119, 7088–7096 (2020)
  * Muehlebach and Jordan [2021] Muehlebach, M., Jordan, M.I.: Optimization with momentum: Dynamical, control-theoretic, and symplectic perspectives. Journal of Machine Learning Research 22(73), 1–50 (2021)
  * Moreau [1988] Moreau, J.J.: Unilateral contact and dry friction in finite freedom dynamics. In: Nonsmooth Mechanics and Applications, pp. 1–82. Springer, Wien (1988)
  * Glocker [2001] Glocker, C.: Set-Valued Force Laws. Springer, Berlin (2001)
  * Studer [2009] Studer, C.: Numerics of Unilateral Contacts and Friction. Springer, Berlin (2009)
  * Wang and Liu [2006] Wang, C., Liu, Q.: Convergence properties of inexact projected gradient methods. Optimization 55(3), 301–310 (2006)
  * Birgin et al. [2003] Birgin, E.G., Martínez, J.M., Raydan, M.: Inexact spectral projected gradient methods on convex sets. IMA Journal of Numerical Analysis 23, 539–559 (2003)
  * Beck and Teboulle [2011] Beck, A., Teboulle, M.: Gradient-based algorithms with applications to signal-recovery problems. In: Convex Optimization in Signal Processing and Communications, pp. 42–88. Cambridge University Press, Cambridge (2011)
  * Bloom et al. [2016] Bloom, V., Griva, I., Quijada, F.: Fast projected gradient method for support vector machines. Optimization and Engineering 17(4), 651–662 (2016)
  * Jaggi [2013] Jaggi, M.: Revisiting Frank-Wolfe: Projection-free sparse convex optimization. Proceedings of Machine Learning Research 28(1), 427–435 (2013)
  * Clarkson [2010] Clarkson, K.L.: Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm. ACM Transactions on Algorithms 6(4), 1–30 (2010)
  * Hazan and Kale [2012] Hazan, E., Kale, S.: Projection-free online learning. Proceeding of the International Conference on Machine Learning, 1–8 (2012)
  * Zhang et al. [2020] Zhang, M., Shen, Z., Mokhtari, A., Hassani, H., Karbasi, A.: One sample stochastic Frank-Wolfe. Proceedings of Machine Learning Research 108, 4012–4023 (2020)
  * Garber and Hazan [2015] Garber, D., Hazan, E.: Faster rates for the Frank-Wolfe method over strongly-convex sets. Proceedings of Machine Learning Research 37, 541–549 (2015)
  * Combettes and Pokutta [2020] Combettes, C.W., Pokutta, S.: Boosting Frank-Wolfe by chasing gradients. Proceedings of Machine Learning Research 119, 2111–2121 (2020)
  * Nemirovski and Yudin [1983] Nemirovski, A.S., Yudin, D.B.: Problem Complexity and Method Efficiency in Optimization. John Wiley & Sons, Chichester (1983)
  * Beck and Teboulle [2003] Beck, A., Teboulle, M.: Mirror descent and nonlinear projected subgradient methods for convex optimization. Operations Research Letters 31(3), 167–175 (2003)
  * Bubeck and Cesa-Bianchi [2012] Bubeck, S., Cesa-Bianchi, N.: Regret analysis of stochastic and nonstochastic multi-armed bandit problems. Foundation and Trends in Machine Learning 5(1), 1–122 (2012)
  * Muehlebach and Jordan [2022] Muehlebach, M., Jordan, M.I.: On constraints in first-order optimization: A view from non-smooth dynamical systems. Journal of Machine Learning Research 23(256), 1–47 (2022)
  * Kolev et al. [2023] Kolev, P., Martius, G., Muehlebach, M.: Online learning under adversarial nonlinear constraints. Advances in Neural Information Processing Systems 36, 1–26 (2023)
  * Schechtman et al. [2023] Schechtman, S., Tiapkin, D., Muehlebach, M., Moulines, É.: Orthogonal directions constrained gradient method: from non-linear equality constraints to Stiefel manifold. Proceedings of Machine Learning Research 195, 1228–1258 (2023)
  * Gonçalves et al. [2017] Gonçalves, D.S., Mucherino, A., Lavor, C., Liberti, L.: Recent advances on the interval distance geometry problem. Journal of Global Optimization 69(3), 525–545 (2017)
  * Ibrahim et al. [2023] Ibrahim, A.A., Muehlebach, M., Bacco, C.D.: Optimal transport with constraints: from mirror descent to classical mechanics. arXiv:2309.04727, 1–14 (2023)
  * Bolte et al. [2018] Bolte, J., Hochart, A., Pauwels, E.: Qualification conditions in semialgebraic programming. SIAM Journal on Optimization 28(2), 1867–1891 (2018)
  * Polyak [1964] Polyak, B.T.: Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics 4(5), 1–17 (1964)
  * Leine and van de Wouw [2008] Leine, R.I., Wouw, N.: Stability and Convergence of Mechanical Systems with Unilateral Constraints. Springer, Berlin (2008)
  * Piazza et al. [2021] Piazza, L.D., Marraffa, V., Satco, B.: Measure differential inclusions: Existence results and minimum problems. Set-Valued and Variational Analysis 29(2), 361–382 (2021)
  * Leine and van de Wouw [2008] Leine, R.I., Wouw, N.: Uniform convergence of monotone measure differential inclusions: With application to the control of mechanical systems with unilateral constraints. International Journal of Bifurcation and Chaos 18(5), 1435–1457 (2008)
  * Nesterov [2004] Nesterov, Y.: Introductory Lectures on Convex Optimization - A Basic Course. Springer, New York (2004)
  * Hastie et al. [2009] Hastie, T., Tibshirani, R., Friedman, J.: The Elements of Statistical Learning, 2nd edn. Springer, New York (2009)
  * Beck and Teboulle [2009] Beck, A., Teboulle, M.: A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences 2(1), 183–202 (2009)
  * Teel [1999] Teel, A.R.: Asymptotic convergence from $\mathcal{L}_{p}$ stability. IEEE Transactions on Automatic Control 44(11), 2169–2170 (1999)