11institutetext: Francisco J. Aragón-Artacho 22institutetext: Department of
Mathematics, University of Alicante, Alicante, Spain
22email<EMAIL_ADDRESS>33institutetext: Boris S. Mordukhovich
44institutetext: Department of Mathematics, Wayne State University, Detroit,
Michigan 48202, USA
44email<EMAIL_ADDRESS>55institutetext: Pedro Pérez-Aros
66institutetext: Instituto de Ciencias de la Ingeniería, Universidad de
O’Higgins, Rancagua, Chile
66email<EMAIL_ADDRESS>
# Coderivative-Based Semi-Newton Method
in Nonsmooth Difference Programming ††thanks: Research of the first author was
partially supported by the Ministry of Science, Innovation and Universities of
Spain and the European Regional Development Fund (ERDF) of the European
Commission, Grant PGC2018-097960-B-C22, and by the Generalitat Valenciana,
grant AICO/2021/165. Research of the second author was partially supported by
the USA National Science Foundation under grants DMS-1808978 and DMS-2204519,
by the Australian Research Council under Discovery Project DP-190100555, and
by the Project 111 of China under grant D21024. Research of the third author
was partially supported by grants: Fondecyt Regular 1190110 and Fondecyt
Regular 1200283.

Francisco J. Aragón-Artacho  Boris S. Mordukhovich Pedro Pérez-Aros

###### Abstract

This paper addresses the study of a new class of nonsmooth optimization
problems, where the objective is represented as a difference of two generally
nonconvex functions. We propose and develop a novel Newton-type algorithm to
solving such problems, which is based on the coderivative generated second-
order subdifferential (generalized Hessian) and employs advanced tools of
variational analysis. Well-posedness properties of the proposed algorithm are
derived under fairly general requirements, while constructive convergence
rates are established by using additional assumptions including the
Kurdyka–Łojasiewicz condition. We provide applications of the main algorithm
to solving a general class of nonsmooth nonconvex problems of structured
optimization that encompasses, in particular, optimization problems with
explicit constraints. Finally, applications and numerical experiments are
given for solving practical problems that arise in biochemical models,
constrained quadratic programming, etc., where advantages of our algorithms
are demonstrated in comparison with some known techniques and results.

###### Keywords:

Nonsmooth difference programming generalized Newton methods global convergence
convergence rates variational analysis generalized differentiation

###### MSC:

49J53, 90C15, 9J52

##  1 Introduction

The primary mathematical model considered in this paper is described by

$\min_{x\in\mathbb{R}^{n}}\varphi(x):=g(x)-h(x),$ (1)

where $g:\mathbb{R}^{n}\to\mathbb{R}$ is of class $\mathcal{C}^{1,1}$ (i.e.,
the collection of $\mathcal{C}^{1}$-smooth functions with locally Lipschitzian
derivatives), and where $h:\mathbb{R}^{n}\to\mathbb{R}$ is a locally
Lipschitzian and prox-regular function; see below. Although (1) is a problem
of unconstrained optimization, it will be shown below that a large class of
constrained optimization problems can be reduced to this form. In what
follows, we label the optimization class in (1) as problems of difference
programming.

The difference form (1) reminds us of problems of DC $($difference of
convex$)$ programming, which have been intensively studied in optimization
with a variety of practical applications; see, e.g., Aragon2020 ; Artacho2019
; AragonArtacho2018 ; Oliveira_2020 ; hiriart ; Toh ; Tao1997 ; Tao1998 ;
Tao1986 and the references therein. However, we are not familiar with a
systematic study of the class of difference programming problems considered in
this paper.

Our main goal here is to develop an efficient numerical algorithm to solve the
class of difference programs (1) with subsequent applications to nonsmooth and
nonconvex problems of particular structures, problems with geometric
constraints, etc. Furthermore, the efficiency of the proposed algorithm and
its modifications is demonstrated by solving some practical models for which
we conduct numerical experiments and compare the obtained results with
previously known developments and computations by using other algorithms. The
proposed algorithm is of a regularized damped Newton type with a novel choice
of directions in the iterative scheme providing a global convergence of
iterates to a stationary point of the cost function. At the first order, the
novelty of our algorithm, in comparison with, e.g., the most popular DCA
algorithm by Tao et al. Tao1997 ; Tao1998 ; Tao1986 and its boosted
developments by Aragón-Artacho et al.Aragon2020 ; Artacho2019 ;
AragonArtacho2018 ; MR4078808 in DC programming, is that instead of a convex
subgradient of $h$ in (1), we now use a limiting subgradient of $-h$. No
second-order information on $h$ is used in what follows. Concerning the other
function $g$ in (1), which is nonsmooth of the second-order, our algorithm
replaces the classical Hessian matrix by the generalized Hessian/second-order
subdifferential of $g$ in the sense of Mordukhovich m92 . The latter
construction, which is defined as the coderivative of the limiting
subdifferential has been well recognized in variational analysis and
optimization due its comprehensive calculus and explicit evaluations for broad
classes of extended-real-valued functions arising in applications. We refer
the reader to, e.g., chhm ; dsy ; Helmut ; hmn ; hos ; hr ;
2020arXiv200910551D ; MR3823783 ; MR2191744 ; mr ; os ; yy and the
bibliographies therein for more details. Note also that the aforementioned
generalized Hessian has already been used in differently designed algorithms
of the Newton type to solve optimization-related problems of different
nonsmooth structures in comparison with (1); see Helmut ; 2020arXiv200910551D
; jogo ; 2021arXiv210902093D ; BorisEbrahim . Having in mind the discussions
above, we label the main algorithm developed in this paper as the regularized
coderivative-based damped semi-Newton method (abbr. RCSN).

The rest of the paper is organized as follows. Section 2 recalls constructions
and statements from variational analysis and generalized differentiation,
which are broadly used in the formulations and proofs of the major results.
Besides well-known facts, we present here some new notions and further
elaborations.

In Section 3, we design our main RCSN algorithm, discuss each of its steps,
and establish various results on its performance depending on imposed
assumptions whose role and importance are illustrated by examples.
Furthermore, Section 4 employs the Kurdyka-Łojasiewicz (KL) property of the
cost function to establish quantitative convergence rates of the RCSN
algorithm depending on the exponent in the KL inequality.

Section 5 addresses the class of (nonconvex) problems of structured
optimization with the cost functions given in the form $f(x)+\psi(x)$, where
$f\colon\mathbb{R}^{n}\to\mathbb{R}$ is a twice continuously differentiable
function with a Lipschitzian Hessian (i.e., of class ${\cal C}^{2,1}$), while
$\psi\colon\mathbb{R}^{n}\to\overline{\mathbb{R}}:=(-\infty,\infty]$ is an
extended-real-valued prox-bounded function. By using the forward-backward
envelope MR3845278 and the associated Asplund function asplund , we reduce
this class of structured optimization problems to the difference form (1) and
then employ the machinery of RCSN to solving problems of this type. As a
particular case of RCSN, we design and justify here a new projected-like
Newton algorithm to solve optimization problems with geometric constraints
given by general closed sets.

Section 6 is devoted to implementations of the designed algorithms and
numerical experiments in two different problems arising in practical modeling.
Although these problems can be treated after some transformations by DCA-like
algorithms, we demonstrate in this section numerical advantages of the newly
designed algorithms over the known developments in both smooth and nonsmooth
settings. The concluding Section 7 summarizes the major achievements of the
paper and discusses some directions of our future research.

##  2 Tools of Variational Analysis and Generalized Differentiation

Throughout the entire paper, we deal with finite-dimensional Euclidean spaces
and use the standard notation and terminology of variational analysis and
generalized differentiation; see, e.g., MR3823783 ; MR1491362 , where the
reader can find the majority of the results presented in this section. Recall
that $\mathbb{B}_{r}(x)$ stands for the closed ball centered at
$x\in\mathbb{R}^{n}$ with radius $r>0$ and that
$\mathbb{N}:=\\{1,2,\ldots\\}$.

Given a set-valued mapping $F:\mathbb{R}^{n}\;{\lower
1.0pt\hbox{$\rightarrow$}}\kern-10.0pt\hbox{\raise
2.0pt\hbox{$\rightarrow$}}\;\mathbb{R}^{m}$, its graph is the set
$\operatorname{gph}F:=\big{\\{}(v,w)\in{\mathbb{R}^{n}}\times\mathbb{R}^{m}\;|\;w\in
F(x)\big{\\}}$, while the (Painlevé–Kuratowski) outer limit of $F$ at
$x\in\mathbb{R}^{n}$ is defined by

$\mathop{{\rm Lim}\,{\rm sup}}_{u\to
x}F(u):=\big{\\{}y\in\mathbb{R}^{m}\;\big{|}\;\exists\,u_{k}\to x,\,y_{k}\to
y,\;y_{k}\in F(u_{k})\;\mbox{as}\;k\in\mathbb{N}\big{\\}}.$ (2)

For a nonempty set $C\subseteq\mathbb{R}^{n}$, the (Fréchet) regular normal
cone and (Mordukhovich) basic/limiting normal cone at $x\in C$ are defined,
respectively, by

$\displaystyle\widehat{N}(x;C)=\widehat{N}_{C}(x):$
$\displaystyle=\Big{\\{}x^{*}\in\mathbb{R}^{n}\;\Big{|}\;\limsup\limits_{u\overset{C}{\to}x}\Big{\langle}x^{\ast},\frac{u-x}{\|u-x\|}\Big{\rangle}\leq
0\Big{\\}},$ (3) $\displaystyle N(x;C)=N_{C}(x):$ $\displaystyle=\mathop{{\rm
Lim}\,{\rm sup}}\limits_{u\overset{C}{\to}x}\widehat{N}(u;C),$

where “$u\overset{C}{\to}x$” means that $u\to x$ with $u\in C$. We use the
convention $\widehat{N}(x;C)=N(x;C):=\emptyset$ if $x\notin C$. The indicator
function $\delta_{C}(x)$ of $C$ is equal to 0 if $x\in C$ and to $\infty$
otherwise.

For a lower semicontinuous (l.s.c.) function
$f:{\mathbb{R}^{n}}\to\overline{\mathbb{R}}$, its domain and epigraph are
given by $\mbox{\rm dom}\,f:=\\{x\in\mathbb{R}^{n}\mid f(x)<\infty\\}$ and
$\mbox{\rm
epi}\,f:=\\{(x,\alpha)\in{\mathbb{R}^{n}}\times\mathbb{R}\;|\;f(x)\leq\alpha\\},$
respectively. The regular and basic subdifferentials of $f$ at $x\in\mbox{\rm
dom}\,f$ are defined by

$\displaystyle\widehat{\partial}f(x)$
$\displaystyle:=\big{\\{}x^{\ast}\in\mathbb{R}^{n}\mid(x^{\ast},-1)\in\widehat{N}\big{(}(x,f(x));\mbox{\rm
epi}\,f\big{)}\big{\\}},$ (4) $\displaystyle\partial f(x)$
$\displaystyle:=\big{\\{}x^{\ast}\in\mathbb{R}^{n}\mid(x^{\ast},-1)\in
N\big{(}(x,f(x));\mbox{\rm epi}\,f\big{)}\big{\\}},$

via the corresponding normal cones (3) to the epigraph. The function $f$ is
said to be lower/subdifferentially regular at $\bar{x}\in\mbox{\rm dom}\,f$ if
$\partial f(\bar{x})=\widehat{\partial}f(\bar{x})$.

Given further a set-valued mapping/multifunction $F:{\mathbb{R}^{n}}\;{\lower
1.0pt\hbox{$\rightarrow$}}\kern-10.0pt\hbox{\raise
2.0pt\hbox{$\rightarrow$}}\;\mathbb{R}^{m}$, the regular and basic
coderivatives of $F$ at $(x,y)\in\operatorname{gph}F$ are defined for all
$y^{*}\in\mathbb{R}^{m}$ via the corresponding normal cones (3) to the graph
of $F$, i.e.,

$\displaystyle\widehat{D}^{\ast}F(x,y)(y^{\ast})$
$\displaystyle:=\big{\\{}x^{\ast}\in{\mathbb{R}^{n}}\;\big{|}\;(x^{\ast},-y^{\ast})\in\widehat{N}\big{(}(x,y);\operatorname{gph}F\big{)}\big{\\}},$
(5) $\displaystyle{D}^{\ast}F(x,y)(y^{\ast})$
$\displaystyle:=\big{\\{}x^{\ast}\in{\mathbb{R}^{n}}\;\big{|}\;(x^{\ast},-y^{\ast})\in
N\big{(}(x,y);\operatorname{gph}F\big{)}\big{\\}},$

where $y$ is omitted if $F$ is single-valued at $x$. When $F$ is single-valued
and locally Lipschitzian around $x$, the basic coderivative has the following
representation via the basic subdifferential of the scalarization

$\displaystyle{D}^{\ast}F(x)(y^{\ast})=\partial\langle
y^{\ast},F\rangle(x),\text{ where }\langle y^{\ast},F\rangle(x):=\langle
y^{\ast},F(x)\rangle.$ (6)

Recall that a set-valued mapping $F:\mathbb{R}^{n}\;{\lower
1.0pt\hbox{$\rightarrow$}}\kern-10.0pt\hbox{\raise
2.0pt\hbox{$\rightarrow$}}\;\mathbb{R}^{m}$ is strongly metrically subregular
at $(\bar{x},\bar{y})\in\operatorname{gph}F$ if there exist
$\kappa,\varepsilon>0$ such that

$\displaystyle\|x-\bar{x}\|\leq\kappa\|y-\bar{y}\|\;\text{ for all
}\;(x,y)\in\mathbb{B}_{\varepsilon}(\bar{x},\bar{y})\cap\operatorname{gph}F.$
(7)

It is well-known that this property of $F$ is equivalent to the calmness
property of the inverse mapping $F^{-1}$ at $(\bar{y},\bar{x})$. In what
follows, we use the calmness property of single-valued mappings
$h\colon\mathbb{R}^{n}\to\mathbb{R}^{m}$ at $\bar{x}$ meaning that there exist
positive numbers $\kappa$ and $\varepsilon>0$ such that

$\|h(x)-h(\bar{x})\|\leq\kappa\|x-\bar{x}\|\;\text{ for all
}\;x\in\mathbb{B}_{\varepsilon}(\bar{x}).$ (8)

The infimum of all $\kappa>0$ in (8) is called the _exact calmness bound_ of
$h$ at $\bar{x}$ and is denoted it by $\mbox{\rm clm}\,h(\bar{x})$. On the
other hand, a multifunction $F\colon\mathbb{R}^{n}\;{\lower
1.0pt\hbox{$\rightarrow$}}\kern-10.0pt\hbox{\raise
2.0pt\hbox{$\rightarrow$}}\;\mathbb{R}^{m}$ is strongly metrically regular
around $(\bar{x},\bar{y})\in\operatorname{gph}F$ if its inverse $F^{-1}$
admits a single-valued and Lipschitz continuous localization around this
point.

Along with the (first-order) basic subdifferential in (4), we consider the
second-order subdifferential/generalized Hessian of
$f:{\mathbb{R}^{n}}\to\overline{\mathbb{R}}$ at $x\in\mbox{\rm dom}\,f$
relative to $x^{\ast}\in\partial f(x)$ defined by

$\partial^{2}f(x,x^{\ast})(v^{\ast})=\left(D^{\ast}\partial
f\right)(x,x^{\ast})(v^{\ast}),\quad v^{\ast}\in{\mathbb{R}^{n}}$ (9)

and denoted by $\partial^{2}f(x)(v^{\ast})$ when $\partial f(x)$ is a
singleton. If $f$ is twice continuously differentiable
($\mathcal{C}^{2}$-smooth) around $x$, then
$\partial^{2}f(x)(v^{\ast})=\\{\nabla^{2}f(x)v^{\ast}\\}$.

Next we introduce an extension of the notion of positive-definiteness for
multifunctions, where the the corresponding constant may not be positive.

###### Definition 1

Let $F:\mathbb{R}^{n}\;{\lower
1.0pt\hbox{$\rightarrow$}}\kern-10.0pt\hbox{\raise
2.0pt\hbox{$\rightarrow$}}\;\mathbb{R}^{n}$ and $\xi\in\mathbb{R}$. Then $F$
is _$\xi$ -lower-definite_ if

$\displaystyle\langle y,x\rangle\geq\xi\|x\|^{2}\;\text{ for all
}\;(x,y)\in\operatorname{gph}F.$ (10)

###### Remark 1

We can easily check the following:

(i) For any symmetric matrix $Q$ with the smallest eigenvalue
$\lambda_{\min}(Q)$, the function $f(x)=Qx$ is $\lambda_{\min}(Q)$-lower-
definite.

(ii) If a function $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is strongly
convex with modulus $\rho>0$, (i.e., $f-\frac{\rho}{2}\|\cdot\|^{2}$ is
convex), it follows from (MR3823783, , Corollary 5.9) that
$\partial^{2}f(x,x^{\ast})$ is $\rho$-lower-definite for all
$(x,x^{\ast})\in\operatorname{gph}\partial f$.

(iii) If $F_{1},F_{2}:\mathbb{R}^{n}\;{\lower
1.0pt\hbox{$\rightarrow$}}\kern-10.0pt\hbox{\raise
2.0pt\hbox{$\rightarrow$}}\;\mathbb{R}^{n}$ are $\xi_{1}$ and $\xi_{2}$-lower-
definite, then the sum $F_{1}+F_{2}$ is $(\xi_{1}+\xi_{2})$-lower-definite.

Recall next that a function $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is
_prox-regular_ at $\bar{x}\in\mathbb{R}^{n}$ _for_ $\bar{v}\in\partial
f(\bar{x})$ if it is l.s.c. around $\bar{x}$ and there exist $\varepsilon>0$
and $r\geq 0$ such that

$\displaystyle f(x^{\prime})\geq f(x)+\langle
v,x^{\prime}-x\rangle-\frac{r}{2}\|x^{\prime}-x\|^{2}$ (11)

whenever $x,x^{\prime}\in\mathbb{B}_{\varepsilon}(\bar{x})$ with $f(x)\leq
f(\bar{x})+\varepsilon$ and $v\in\partial
f(x)\cap\mathbb{B}_{\varepsilon}(\bar{v})$. If this holds for all
$\bar{v}\in\partial f(\bar{x})$, $f$ is said to be _prox-regular at_
$\bar{x}$.

###### Remark 2

The class of prox-regular functions has been well-recognized in modern
variational analysis. It is worth mentioning that if $f$ is a locally
Lipschitzian function around $\bar{x}$, then the following properties of $f$
are equivalent: (i) prox-regularity at $\bar{x}$, (ii) lower-$\mathcal{C}^{2}$
at $\bar{x}$, and (iii) primal-lower-nice at $\bar{x}$; see, e.g., (MR2101873,
, Corollary 3.12) for more details.

Given a function $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ and
$\bar{x}\in\mbox{\rm dom}\,f$, the upper directional derivative of $f$ at
$\bar{x}$ with respect to $d\in\mathbb{R}^{n}$ is defined by

$f^{\prime}(\bar{x};d):=\limsup\limits_{t\to
0^{+}}\frac{f(\bar{x}+td)-f(\bar{x})}{t}.$ (12)

The following proposition establishes various properties of prox-regular
functions used below. We denote the convex hull of a set by “co”.

###### Proposition 1

Let $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ be locally Lipschitzian around
$\bar{x}$ and prox-regular at this point. Then $f$ is lower regular at
$\bar{x}$, $\mbox{\rm co}\,\partial(-f)(\bar{x})=-\partial f(\bar{x})$, and
for any $d\in\mathbb{R}^{n}$ we have the representations

$\displaystyle(-f)^{\prime}(\bar{x};d)=\inf\big{\\{}\langle
w,d\rangle\;\big{|}\;w\in\partial(-f)(\bar{x})\big{\\}}=\inf\big{\\{}\langle
w,d\rangle\;\big{|}\;w\in-\partial f(\bar{x})\big{\\}}.$ (13)

###### Proof

First we fix an arbitrary subgradient $\bar{v}\in\partial f(\bar{x})$ and
deduce from (11) applied to $x=\bar{x}$ and $v=\bar{v}$ that

$f(x^{\prime})\geq
f(\bar{x})+\langle\bar{v},x^{\prime}-\bar{x}\rangle-\frac{r}{2}\|x^{\prime}-\bar{x}\|^{2}\;\text{
for all }\;x^{\prime}\in\mathbb{B}_{\varepsilon}(\bar{x}).$

Passing to the limit as $x^{\prime}\to\bar{x}$ tells us that

$\displaystyle\liminf_{x^{\prime}\to\bar{x}}\frac{f(x^{\prime})-f(\bar{x})-\langle\bar{v},x^{\prime}-\bar{x}\rangle}{\|x^{\prime}-\bar{x}\|}\geq
0,$

which means that $\bar{v}\in\widehat{\partial}f(\bar{x})$ and thus shows that
$f$ is lower regular at $\bar{x}$. By the Lipschitz continuity of $f$ around
$\bar{x}$ and the convexity of the set $\widehat{\partial}f(\bar{x})$, we have
that $\widehat{\partial}f(\bar{x})=\partial f(\bar{x})=\mbox{\rm co}\,\partial
f(\bar{x})=\overline{\partial}f(\bar{x})$, where $\overline{\partial}$ denotes
the (Clarke) generalized gradient. It follows from
$\overline{\partial}(-f)(\bar{x})=-\overline{\partial}f(\bar{x})$ that
$\overline{\partial}(-f)(\bar{x})=-\partial f(\bar{x})$, which implies
therefore that $\partial(-f)(\bar{x})\subseteq-\partial f(\bar{x})$.

Pick $v\in\partial(-f)(\bar{x})$, $d\in\mathbb{R}^{n}$ and find by the prox-
regularity of $f$ at $\bar{x}$ for $-v\in\partial f(\bar{x})$ that there
exists $r>0$ such that

$\displaystyle\langle
v,d\rangle+\frac{rt}{2}\|d\|^{2}\geq\frac{-f(\bar{x}+td)+f(\bar{x})}{t}$

if $t>0$ is small enough. This yields $(-f)^{\prime}(\bar{x};d)\leq\langle
v,d\rangle$ for all ${v}\in\partial(-f)(\bar{x})$ and thus verifies the
inequality “$\leq$” in the first representation of (13).

To prove the opposite inequality therein, take $t_{k}\to 0^{+}$ such that

$\displaystyle\lim_{k\to\infty}\frac{-f(\bar{x}+t_{k}d)+f(\bar{x})}{t_{k}}=(-f)^{\prime}(\bar{x};d).$

Employing the mean value theorem from (MR3823783, , Corollary 4.12)) gives us

$\displaystyle f(\bar{x}+t_{k}d)-f(\bar{x})=t_{k}\langle v_{k},d\rangle\text{
for some }v_{k}\in\partial f(\bar{x}+\lambda_{k}t_{k}d)\text{ with
}\lambda_{k}\in(0,1).$

It follows from the Lipschitz continuity of $f$ that $\\{v_{k}\\}$ is bounded,
and so we can assume that $v_{k}\to\bar{v}\in\partial f(\bar{x})$. Therefore,

$\begin{array}[]{ll}(-f)^{\prime}(\bar{x};d)=\langle-\bar{v},d\rangle\geq\inf\big{\\{}\langle
w,d\rangle\;\big{|}\;w\in-\partial f(\bar{x})\big{\\}}\\\
=\inf\big{\\{}\langle w,d\rangle\;\big{|}\;w\in\mbox{\rm
co}\,\partial(-f)(\bar{x})\big{\\}}=\inf\big{\\{}\langle
w,d\rangle\;\big{|}\;w\in\partial(-f)(\bar{x})\big{\\}},\end{array}$

which verifies (13) and completes the proof of the proposition.

Next we define the notion of stationary points for problem (1) the finding of
which is the goal of our algorithms.

###### Definition 2

Let $\varphi=g-h$ be the cost function in (1), where $g$ is of class
$\mathcal{C}^{1,1}$ around some point $\bar{x}$, and where $h$ is locally
Lipschitzian around $\bar{x}$ and prox-regular at this point. Then $\bar{x}$
is a _stationary point_ of (1) if $0\in\partial\varphi(\bar{x})$.

###### Remark 3

The stationarity notion $0\in\partial\varphi(\bar{x})$, expressed via the
limiting subdiffential, is known as the M$($ordukhovich$)$-stationarity. Since
no other stationary points are considered in this paper, we skip “M” in what
follows. Observe from the subdifferential sum rule in our setting that
$\bar{x}$ is a stationary point in (1) if and only if $0\in\nabla
g(\bar{x})+\partial(-h)(\bar{x})$. Thus every stationary point $\bar{x}$ is a
critical point in the sense that $0\in\nabla g(\bar{x})-\partial h(\bar{x})$.
By Proposition 1, the latter can be equivalently described in terms of the
generalized gradient and also via the symmetric subdifferential MR3823783 of
$\varphi$ at $\bar{x}$ defined by

$\partial^{0}\varphi(\bar{x}):=\partial\varphi(\bar{x})\cup\big{(}-\partial(-\varphi)(\bar{x})\big{)}$
(14)

which possesses the plus-minus symmetry
$\partial^{0}(-\varphi(\bar{x}))=-\partial^{0}(\varphi(\bar{x}))$. When both
$g$ and $h$ are convex, the classical DC algorithm Tao1986 ; Tao1997 and its
BDCA variant MR4078808 can be applied for solving problem (1). Although these
algorithms only converge to critical points, they can be easily combined as in
Aragon2020 with a basic derivative-free optimization scheme to converge to
d-stationary points, which satisfy $\partial h(\bar{x})=\\{\nabla
g(\bar{x})\\}$ (or, equivalently, $\varphi^{\prime}(\bar{x};d)=0$ for all
$d\in\mathbb{R}^{n}$; see (Aragon2020, , Proposition 1)). In the DC setting,
every local minimizer of problem (1) is a d-stationary point (Toland1979, ,
Theorem 3), a property which is stronger than the notion of stationarity in
Definition 2.

To proceed, recall that a mapping $f:U\to\mathbb{R}^{m}$ defined on an open
set $U\subseteq\mathbb{R}^{n}$ is _semismooth_ at $\bar{x}$ if it is locally
Lipschitzian around $\bar{x}$, directionally differentiable at this point, and
the limit

$\displaystyle\lim\limits_{A\in\tiny{\mbox{\rm
co}\,}\overline{\nabla}f(\bar{x}+tu^{\prime}),\atop u^{\prime}\to u,t\to
0^{+}}Au^{\prime}$

exists for all $u\in\mathbb{R}^{n}$, where
$\overline{\nabla}f(x):=\\{A\;|\;\exists x_{k}\overset{D}{\to}x\text{ and
}\nabla f(x_{k})\to A\\}$, and where $D$ is the set on which $f$ is
differentiable; see MR1955649 ; MR3289054 for more details. We say that a
function $g:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is semismoothly
differentiable at $\bar{x}$ if $g$ is $\mathcal{C}^{1}$-smooth around
$\bar{x}$ and its gradient mapping $\nabla g$ is semismooth at this point.

Recall further that a function $\psi:\mathbb{R}^{n}\to\overline{\mathbb{R}}$
is _prox-bounded_ if there exists $\lambda>0$ such that
${\mathtt{e}}_{\lambda}\psi(x)>-\infty$ for some $x\in\mathbb{R}^{n}$, where
${\mathtt{e}}_{\lambda}\psi:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is the
_Moreau envelope_ of $\psi$ with parameter $\lambda>0$ defined by

${\mathtt{e}}_{\lambda}\psi(x):=\inf_{z\in\mathbb{R}^{n}}\Big{\\{}\psi(z)+\frac{1}{2\lambda}\|x-z\|^{2}\Big{\\}}.$
(15)

The number
$\lambda_{\psi}:=\sup\\{\lambda>0\;|\;{\mathtt{e}}_{\lambda}\psi(x)>-\infty\text{
for some }x\in\mathbb{R}^{n}\\}$ is called the _threshold_ of prox-boundedness
of $\psi$. The corresponding _proximal mapping_ is the multifunction
${\mathtt{Prox}}_{\lambda\psi}:\mathbb{R}^{n}\;{\lower
1.0pt\hbox{$\rightarrow$}}\kern-10.0pt\hbox{\raise
2.0pt\hbox{$\rightarrow$}}\;\mathbb{R}^{n}$ given by

${\mathtt{Prox}}_{\lambda\psi}(x):=\mathop{\rm
argmin}_{z\in\mathbb{R}^{n}}\Big{\\{}\psi(z)+\frac{1}{2\lambda}\|x-z\|^{2}\Big{\\}}.$
(16)

Next we observe that that the Moreau envelope can be represented as a DC
function. For any function $\varphi:\mathbb{R}^{n}\to\overline{\mathbb{R}}$,
consider its _Fenchel conjugate_

$\phi^{*}(x):=\sup_{z\in\mathbb{R}^{n}}\big{\\{}\langle
x,z\rangle-\phi(z)\big{\\}},$

and for any $\psi\colon\mathbb{R}^{n}\to\overline{\mathbb{R}}$ and
$\lambda>0$, define the Asplund function

${{\mathtt{A}}_{\lambda}\psi}(x):=\sup\limits_{z\in\mathbb{R}^{n}}\Big{\\{}\frac{1}{\lambda}\langle
z,x\rangle-\psi(z)-\frac{1}{2\lambda}\|z\|^{2}\Big{\\}}=\Big{(}\psi+\frac{1}{2\lambda}\|\cdot\|^{2}\Big{)}^{\ast}(x),$
(17)

which is inspired by Asplund’s study of metric projections in asplund . The
following proposition presents the precise formulation of the aforementioned
statement and reveals some remarkable properties of the Asplund function (17).

###### Proposition 2

Let $\psi$ be a prox-bounded function with threshold $\lambda_{\psi}$. Then
for every $\lambda\in(0,\lambda_{\psi})$, we have the representation

${\mathtt{e}}_{\lambda}\psi(x)=\frac{1}{2\lambda}\|x\|^{2}-{{\mathtt{A}}_{\lambda}\psi}(x),\quad
x\in\mathbb{R}^{n},$ (18)

where the Asplund function is convex and Lipschitz continuous on
$\mathbb{R}^{n}$. Furthermore, for any $x\in\mathbb{R}^{n}$ the following
subdifferential evaluations hold:

$\displaystyle{}\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)$
$\displaystyle\subseteq-\frac{1}{\lambda}{\mathtt{Prox}}_{\lambda\psi}(x),$
(19) $\displaystyle\partial{{\mathtt{A}}_{\lambda}\psi}(x)$
$\displaystyle=\frac{1}{\lambda}\mbox{\rm
co}\,\left({\mathtt{Prox}}_{\lambda\psi}(x)\right).$ (20)

Moreover, if $v\in{\mathtt{Prox}}_{\lambda\psi}(x)$ is such that
$v\notin\mbox{\rm
co}\,\left({\mathtt{Prox}}_{\lambda\psi}(x)\backslash\\{v\\}\right)$, then the
vector $-\frac{1}{\lambda}v$ belongs to
$\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)$. If in addition $f$ is of class
$\mathcal{C}^{2,1}$ on $\mathbb{R}^{n}$, then the function
$x\mapsto{{\mathtt{A}}_{\lambda}\psi}(x-\lambda\nabla f(x))$ is prox-regular
at any point $x\in\mathbb{R}^{n}$.

###### Proof

Representation (18) easily follows from definitions of the Moreau envelope and
Asplund function. Due to the second equality in (17), the Asplund function is
convex on $\mathbb{R}^{n}$. It is also Lipschitz continuous due its finite-
valuedness on $\mathbb{R}^{n}$, which is induced by this property of the
Moreau envelope. The subdifferential evaluations in (19) and (MR1491362, ,
Example 10.32) and the subdifferential sum rule in (MR3823783, , Proposition
1.30)) tell us that
$\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)=-\lambda^{-1}x+\partial({\mathtt{e}}_{\lambda}\psi)(x)$
and
$\partial{{\mathtt{A}}_{\lambda}\psi}(x)=\lambda^{-1}x+\partial(-{\mathtt{e}}_{\lambda}\psi)(x)$
for any $x\in\mathbb{R}^{n}$.

Take further $v\in{\mathtt{Prox}}_{\lambda\psi}(x)$ with
$-\frac{1}{\lambda}v\not\in\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)$ and
show that $v\in{\rm
co}\left({\mathtt{Prox}}_{\lambda\psi}(x)\backslash\\{v\\}\right)$. Indeed, it
follows from (19) that

$\displaystyle\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)$
$\displaystyle\subseteq-\frac{1}{\lambda}{\mathtt{Prox}}_{\lambda\psi}(x)\backslash\\{v\\}.$

The Lipschitz continuity and convexity of ${{\mathtt{A}}_{\lambda}\psi}$
implies that

$\displaystyle{\rm
co}\,\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)=-\partial{{\mathtt{A}}_{\lambda}\psi}(x)$
(21)

by (MR2191744, , Theorem 3.57), which allows us to deduce from (20) and (21)
that

${\rm co}\big{(}{\mathtt{Prox}}_{\lambda\psi}(x)\big{)}={\rm
co}\big{(}{\mathtt{Prox}}_{\lambda\psi}(x)\backslash\\{v\\}\big{)}.$

This verifies the inclusion $v\in\mbox{\rm
co}\,\left({\mathtt{Prox}}_{\lambda\psi}(x)\backslash\\{v\\}\right)$ as
claimed.

Observe finally that the function
$x\mapsto{{\mathtt{A}}_{\lambda}\psi}_{\lambda}(x-\lambda\nabla f(x))$ is the
composition of the convex function ${{\mathtt{A}}_{\lambda}\psi}$ and the
$\mathcal{C}^{1,1}$ mapping $x\mapsto x-\lambda\nabla f(x)$, which ensures by
(MR2069350, , Proposition 2.3) its prox-regularity at any point
$x\in\mathbb{R}^{n}$.

The following remark discusses a useful representation of the basic
subdifferential of the function $-{{\mathtt{A}}_{\lambda}\psi}$ and other
functions of this type.

###### Remark 4

It is worth mentioning that the subdifferential
$\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)$ can be expressed via the set
$D:=\\{x\in\mathbb{R}^{n}\;|\;{{\mathtt{A}}_{\lambda}\psi}\;\mbox{ is
differentiable at }\;x\\}$ as follows:

$\partial(-{{\mathtt{A}}_{\lambda}\psi})(x)=\big{\\{}v\in\mathbb{R}^{n}\;\big{|}\;\text{
there exists }x_{k}\overset{D}{\to}x\text{ and
}\nabla{{\mathtt{A}}_{\lambda}\psi}(x_{k})\to-v\big{\\}}.$ (22)

We refer to (MR1491362, , Theorem 10.31) for more details. Note that we do not
need to take the convex hull on the right-hand side of (22) as in the case of
the generalized gradient of locally Lipschitzian functions.

Finally, recall the definitions of the convergence rates used in the paper.

###### Definition 3

Let $\\{x_{k}\\}$ be a sequence in $\mathbb{R}^{n}$ converging to $\bar{x}$ as
$k\rightarrow\infty$. The convergence rate is said to be:

(i) _R-linear_ if there exist $\mu\in(0,1),c>0$, and $k_{0}\in\mathbb{N}$ such
that

$\left\|x_{k}-\bar{x}\right\|\leq c\mu^{k}\;\text{ for all }\;k\geq k_{0}.$

(ii) _Q-linear_ if there exists $\mu\in(0,1)$ such that

$\limsup_{k\to\infty}\frac{\left\|x_{k+1}-\bar{x}\right\|}{\left\|x_{k}-\bar{x}\right\|}=\mu.$

(iii) _Q-superlinear_ if it is Q-linear for all $\mu\in(0,1)$, i.e., if

$\lim_{k\to\infty}\frac{\left\|x_{k+1}-\bar{x}\right\|}{\left\|x_{k}-\bar{x}\right\|}=0.$

(iv) _Q-quadratic_ if we have

$\limsup_{k\to\infty}\frac{\left\|x_{k+1}-\bar{x}\right\|}{\left\|x_{k}-\bar{x}\right\|^{2}}<\infty.$

##  3 Regularized Coderivative-Based Damped Semi-Newton Method in Nonsmooth
Difference Programming

The goal of this section is to justify the well-posedness and good performance
of the novel algorithm RCSN under appropriate and fairly general assumptions.
In the following remark, we discuss the difference between the choice of
subgradients and hence of directions in RCSN and DC algorithms.

Our main RCSN algorithm to find stationary points of nonsmooth problems (1) of
difference programming is labeled below as Algorithm 1.

1:$x_{0}\in\mathbb{R}^{n}$, $\beta\in(0,1)$, $\zeta>0$, $t_{\min}>0$,
$\rho_{\max}>0$ and $\sigma\in(0,1)$.

2:for $k=0,1,\ldots$ do

3: Take $w_{k}\in\partial\varphi(x_{k})$. If $w_{k}=0$, STOP and return
$x_{k}$.

4: Choose $\rho_{k}\in[0,\rho_{\max}]$ and
$d_{k}\in\mathbb{R}^{n}\backslash\\{0\\}$ such that  $\displaystyle-
w_{k}\in\partial^{2}g(x_{k})(d_{k})+\rho_{k}d_{k}\quad\text{and}\quad\langle
w_{k},d_{k}\rangle\leq-\zeta\|d_{k}\|^{2}.$ (23)

5: Choose any $\overline{\tau}_{k}\geq t_{\min}$. Set
$\overline{\tau}_{k}:=\tau_{k}$.

6: while $\varphi(x_{k}+\tau_{k}d_{k})>\varphi(x_{k})+\sigma\tau_{k}\langle
w_{k},d_{k}\rangle$ do

7: $\tau_{k}=\beta\tau_{k}$.

8: end while

9: Set $x_{k+1}:=x_{k}+\tau_{k}d_{k}$.

10:end for

Algorithm 1 Regularized coderivative-based damped semi-Newton algorithm for
nonsmooth difference programming

###### Remark 5

Observe that Step 2 of Algorithm 1 selects
$w_{k}\in\partial\varphi(x_{k})=\nabla g(x_{k})+\partial(-h)(x_{k})$, which is
equivalent to choosing $v_{k}:=w_{k}-\nabla g(x_{k})$ in the basic
subdifferential of $-h$ at $x_{k}$. Under our assumptions, the set
$\partial(-h)(x_{k})$ can be considerably smaller than $\partial h(x_{k})$;
see the proof of Proposition 1 and also Remark 4 above. Therefore, Step 2
differs from those in DC algorithms, which choose subgradients in $\partial
h(x_{k})$. The purpose of our development is to find a stationary point
instead of a (classical) critical point for problem (1). In some applications,
Algorithm 1 would not be implementable if the user only has access to
subgradients contained in $\partial h(x_{k})$ instead of
$\partial(-h)(x_{k})$. In such cases, a natural alternative to Algorithm 1
would be a scheme replacing $w_{k}\in\partial\varphi(x_{k})$ in Step 2 by
$w_{k}:=\nabla g(x_{k})+v_{k}$ with $v_{k}\in\partial h(x_{k})$. Under the
setting of our convergence results, the modified algorithm would find a
critical point for problem (1), which is not guaranteed to be stationary.

The above discussions are illustrated by the following example.

###### Example 1

Consider problem (1) with $g(x):=\frac{1}{2}x^{2}$ and $h(x):=|x|$. If an
algorithm similar to Algorithm 1 was run by using $x_{0}=0$ as the initial
point but choosing $w_{0}=\nabla g(x_{0})+v_{0}$ with $v_{0}=0\in\partial
h(0)$ (instead of $w_{0}\in\partial\varphi(x_{0})$), it would stop at the
first iteration and return $x=0$, which is a critical point, but not a
stationary one. On the other hand, for any
$w_{0}\in\partial\varphi(0)=\\{-1,1\\}$ we get $w_{0}\neq 0$, and so Algorithm
1 will continue iterating until it converges to one of the two stationary
points $-1/2$ and $1/2$, which is guaranteed by our main convergence result;
see Theorem 3.1 below.

The next lemma shows that Algorithm 1 is well-defined by proving the existence
of a direction $d_{k}$ satisfying (23) in Step 3 for sufficiently large
regularization parameters $\rho_{k}$.

###### Lemma 1

Let $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ be the objective function in problem
(1) with $g\in\mathcal{C}^{1,1}$ and $h$ being locally Lipschitz around
$\bar{x}$ and prox-regular at this point. Further, assume that
$\partial^{2}g(\bar{x})$ is $\xi$-lower-definite for some $\xi\in\mathbb{R}$
and consider a nonzero subgradient $w\in\partial\varphi(\bar{x})$. Then for
any $\zeta>0$ and any $\rho\geq\zeta-\xi$, there exists a nonzero direction
$d\in\mathbb{R}^{n}$ satisfying the inclusion

$\displaystyle-w\in\partial^{2}g(\bar{x})(d)+\rho d.$ (24)

Moreover, any nonzero direction from (24) obeys the conditions:
(i) $\varphi^{\prime}(\bar{x};d)\leq\langle w,d\rangle\leq-\zeta\|d\|^{2}$.
(ii) Whenever $\sigma\in(0,1)$, there exists $\eta>0$ such that

$\displaystyle\varphi(\bar{x}+\tau d)<\varphi(\bar{x})+\sigma\tau\langle
w,d\rangle\leq\varphi(\bar{x})-\sigma\zeta\tau\|d\|^{2}\;\mbox{ when
}\;\tau\in(0,\eta).$

###### Proof

Consider the function $\psi(x):=g(x)+\langle w-\nabla
g(\bar{x}),x\rangle+\frac{\rho}{2}\|x\|^{2}$ for which we clearly have that
$\partial^{2}\psi(\bar{x})=\partial^{2}g(\bar{x})+\rho I$, where $I$ denotes
the identity mapping. This shows by Remark 1 that $\partial^{2}\psi(\bar{x})$
is $(\xi+\rho)$-lower-definite, and thus it is $\zeta$-lower-definite as well.
Since $\nabla\psi(\bar{x})=w\neq 0$ and $\zeta>0$, it follows from
(2021arXiv210902093D, , Proposition 3.1) (which requires $\psi$ to be
$\mathcal{C}^{1,1}$ on $\mathbb{R}^{n}$, but actually only $\mathcal{C}^{1,1}$
around $\bar{x}$ is needed) that there exists a nonzero direction $d$ such
that $-\nabla\psi(\bar{x})\in\partial^{2}\psi(\bar{x})(d)$. This readily
verifies (24), which yields in turn the second inequality in (i) due to
Definition 1. On the other hand, we have by Proposition 1 the following:

$\displaystyle\begin{aligned} \varphi^{\prime}(\bar{x};d)=\lim\limits_{t\to
0^{+}}\frac{g(\bar{x}+td)-g(\bar{x})}{t}+\limsup\limits_{t\to
0^{+}}\frac{-h(\bar{x}+td)+h(\bar{x})}{t}\\\ =\langle\nabla
g(\bar{x}),d\rangle+\inf\big{\\{}\langle w,d\rangle\;\big{|}\;w\in-\partial
h(\bar{x})\big{\\}}\leq\langle\nabla
g(\bar{x})+v,d\rangle\leq-\zeta\|d\|^{2},\end{aligned}$ (25)

where in the last estimate is a consequence of the second inequality in (i).

Finally, assertion (ii) follows directly from (25) and the definition of
directional derivatives (12).

###### Remark 6

Under the $\xi$-lower-definiteness of $\partial^{2}g(x_{k})$, Lemma 1
guarantees the existence of a direction $d_{k}$ satisfying both conditions in
(23) for all $\rho_{k}\geq\zeta-\xi$. When $\xi$ is unknown, it is still
possible to implement Step 3 of the algorithm as follows. Choose first any
initial value of $\rho\geq 0$, then compute a direction satisfying the
inclusion in (23) and continue with Step 4 if the descent condition in (23)
holds. Otherwise, increase the value of $\rho$ and repeat the process until
the descent condition is satisfied.

The next example demonstrates that the prox-regularity of $h$ is not a
superfluous assumption in Lemma 1. Namely, without it the direction $d$ used
in Step 3 of Algorithm 1 can even be an ascent direction.

###### Example 2

Consider the least squares problem given by

$\min_{x\in\mathbb{R}^{2}}\frac{1}{2}(Ax-b)^{2}+\|x\|_{1}-\|x\|_{2},\quad
x\in\mathbb{R}^{2},$

with $A:=[1,0]$ and $b:=1$. Denote $g(x):=\frac{1}{2}\|Ax-b\|^{2}$ and
$h(x):=\|x\|_{2}-\|x\|_{1}$. If we pick $\bar{x}:=(1,0)^{T}$, the function $h$
is not prox-regular at $\bar{x}$ because it is not lower regular at $\bar{x}$;
see Proposition 1. Indeed, $\widehat{\partial}h(\bar{x})=\emptyset$, while

$\partial
h(\bar{x})=\frac{\bar{x}}{\|\bar{x}\|}+\partial(-\|\cdot\|_{1})(\bar{x})=\left\\{\begin{pmatrix}0\\\
-1\end{pmatrix},\begin{pmatrix}0\\\ 1\end{pmatrix}\right\\}.$

Therefore, although $\nabla^{2}g(\bar{x})=A^{T}A$ is
$\lambda_{\min}(A^{T}A)$-lower-definite, the assumptions of Lemma 1 are not
satisfied. Due to the representation

$\partial(-h)(\bar{x})=-\frac{\bar{x}}{\|\bar{x}\|}+\partial\|\cdot\|_{1}(\bar{x})=\left\\{\begin{pmatrix}0\\\
v\end{pmatrix}\;\Bigg{|}\;v\in[-1,1]\right\\},$

the choice of $v:=(0,1)^{T}\in\partial(-h)(\bar{x})$ yields $w:=\nabla
g(\bar{x})+v=(0,1)^{T}\in\partial\varphi(\bar{x})$. For any $\rho>0$,
inclusion (24) gives us $d=(0,-1/\rho)^{T}$. This is an ascent direction for
the objective function $\varphi(x)=g(x)-h(x)$ at $\bar{x}$ due to

$\varphi(\bar{x}+\tau
d)=1+\frac{\tau}{\rho}-\sqrt{1+(\tau/\rho)^{2}}>\varphi(\bar{x})=0\;\mbox{ for
all }\;\tau>0,$

which illustrates that the prox-regularity is an essential assumption in Lemma
1.

Algorithm 1 either stops at a stationary point, or produces an infinite
sequence of iterates. The convergence properties of the iterative sequence of
our algorithm are obtained below in the main theorem of this section. Prior to
the theorem, we derive yet another lemma, which establishes the following
descent property for the difference of a $\mathcal{C}^{1,1}$ function and a
prox-regular one.

###### Lemma 2

Let $\varphi(x)=g(x)-h(x)$, where $g$ is of class $\mathcal{C}^{1,1}$ around
$\bar{x}$, and where $h$ is continuous around $\bar{x}$ and prox-regular at
this point. Then for every $\bar{v}\in\partial h(\bar{x})$, there exist
positive numbers $\varepsilon$ and $r$ such that

$\displaystyle\varphi(y)\leq\varphi(x)+\langle\nabla
g(x)-v,y-x\rangle+r\|y-x\|^{2}$

whenever $x,y\in\mathbb{B}_{\varepsilon}(\bar{x})$ and $v\in\partial
h(x)\cap\mathbb{B}_{\varepsilon}(\bar{v})$.

###### Proof

Pick any $\bar{v}\in\partial h(\bar{x})$ and deduce from the imposed prox-
regularity and continuity of $h$ that there exist $\varepsilon_{1}>0$ and
$r_{1}>0$ such that

$\displaystyle-h(y)\leq-h(x)+\langle-v,y-x\rangle+r_{1}\|y-x\|^{2}\;\mbox{ for
all }\;x,y\in\mathbb{B}_{\varepsilon_{1}}(\bar{x})$ (26)

and all $v\in\partial h(x)\cap\mathbb{B}_{\varepsilon_{1}}(\bar{v})$. It
follows from the $\mathcal{C}^{1,1}$ property of $g$ by (MR3289054, , Lemma
A.11) that there exist positive numbers $r_{2}$ and $\varepsilon_{2}$ such
that

$\displaystyle g(y)\leq g(x)+\langle\nabla
g(x),y-x\rangle+r_{2}\|y-x\|^{2}\;\mbox{ for all
}\;\mathbb{B}_{\varepsilon_{2}}.$ (27)

Summing up the inequalities in (26) and (27) and defining $r:=r_{1}+r_{2}$ and
$\varepsilon:=\min\\{\varepsilon_{1},\varepsilon_{2}\\}$, we get that

$\displaystyle g(y)-h(y)\leq g(x)-h(x)+\langle\nabla
g(x)-v,y-x\rangle+r\|y-x\|^{2}$

for all $x,y\in\mathbb{B}_{\varepsilon}(\bar{x})$ and all $v\in\partial
h(x)\cap\mathbb{B}_{\varepsilon}(\bar{v})$. This completes the proof.

Now we are ready to establish the aforementioned theorem about the performance
of Algorithm 1.

###### Theorem 3.1

Let $\varphi:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ be the objective function
of problem (1) given by $\varphi=g-h$ with $\inf\varphi>-\infty$. Pick an
initial point $x_{0}\in\mathbb{R}^{n}$ and suppose that the sublevel set
$\Omega:=\\{x\in\mathbb{R}^{n}\;|\;\varphi(x)\leq\varphi(x_{0})\\}$ is closed.
Assume also that:
(a) The function $g$ is $\mathcal{C}^{1,1}$ around every $x\in\Omega$ and the
second-order subdifferential $\partial^{2}g(x)$ is $\xi$-lower-definite for
all $x\in\Omega$ with some $\xi\in\mathbb{R}$.
(b) The function $h$ is locally Lipschitzian and prox-regular on $\Omega$.
Then Algorithm 1 either stops at a stationary point, or produces sequences
$\\{x_{k}\\}\subseteq\Omega$, $\\{\varphi(x_{k})\\}$, $\\{w_{k}\\}$,
$\\{d_{k}\\}$, and $\\{\tau_{k}\\}$ such that:
(i) The sequence $\\{\varphi(x_{k})\\}$ monotonically decreases and converges.
(ii) If $\\{x_{k_{j}}\\}$ as $j\in\mathbb{N}$ is any bounded subsequence of
$\\{x_{k}\\}$, then $\displaystyle\inf_{j\in\mathbb{N}}\tau_{k_{j}}>0$,

$\displaystyle\sum\limits_{j\in\mathbb{N}}\|d_{k_{j}}\|^{2}<\infty,\;\sum\limits_{j\in\mathbb{N}}\|x_{k_{j}+1}-x_{k_{j}}\|^{2}<\infty,\;\text{
and }\;\sum\limits_{j\in\mathbb{N}}\|w_{k_{j}}\|^{2}<\infty.$

In particular, the boundedness of the entire sequence $\\{x_{k}\\}$ ensures
that the set of accumulation points of $\\{x_{k}\\}$ is a nonempty, closed,
and connected.
(iii) If $x_{k_{j}}\to\bar{x}$ as $j\to\infty$, then $\bar{x}$ is a stationary
point of problem (1) with the property
$\varphi(\bar{x})=\displaystyle\inf_{k\in\mathbb{N}}\varphi(x_{k})$.
(iv) If $\\{x_{k}\\}$ has an isolated accumulation point $\bar{x}$, then the
entire sequence $\\{x_{k}\\}$ converges to $\bar{x}$ as $k\to\infty$, where
$\bar{x}$ is a stationary point of (1).

###### Proof

If Algorithm 1 stops after a finite number of iterations, then it clearly
returns a stationary point. Otherwise, it produces an infinite sequence
$\\{x_{k}\\}$. By Step 5 of Algorithm 1 and Lemma 1, we have that
$\inf\varphi\leq\varphi(x_{k+1})<\varphi(x_{k})$ for all $k\in\mathbb{N}$,
which proves assertion (i) and also shows that $\\{x_{k}\\}\subseteq\Omega$.

To proceed, suppose that $\\{x_{k}\\}$ has a bounded subsequence
$\\{x_{k_{j}}\\}$ (otherwise there is nothing to prove) and split the rest of
the proof into the five claims.

Claim 1: _The sequence $\\{\tau_{k_{j}}\\}$, associated with $\\{x_{k_{j}}\\}$
as $j\in\mathbb{N}$ and produced by Algorithm 1, is bounded from below._
Indeed, otherwise consider a subsequence $\\{\tau_{\nu_{i}}\\}$ of
$\\{\tau_{k_{j}}\\}$ such that $\tau_{\nu_{i}}\to 0^{+}$ as $i\to\infty$.
Since $\\{x_{k_{j}}\\}$ is bounded, we can assume that $\\{x_{\nu_{i}}\\}$
converges to some point $\bar{x}$. By Lemma 1, we have that

$\displaystyle-\langle
w_{\nu_{i}},d_{\nu_{i}}\rangle\geq\zeta\|d_{\nu_{i}}\|^{2}\;\mbox{ for all
}\;i\in\mathbb{N},$ (28)

which yields by the Cauchy–Schwarz inequality the estimate

$\displaystyle\|w_{\nu_{i}}\|\geq\zeta\|d_{\nu_{i}}\|,\quad i\in\mathbb{N}.$
(29)

Since $\varphi$ is locally Lipschitzian and
$w_{\nu_{i}}\in\partial\varphi(x_{{\nu_{i}}})$, we suppose without loss of
generality that $w_{\nu_{i}}$ converges to some
$\bar{w}\in\partial\varphi(\bar{x})\subseteq\nabla g(\bar{x})-\partial
h(\bar{x})$ as $i\to\infty$. It follows from (29) that $\\{d_{\nu_{i}}\\}$ is
bounded, and therefore $d_{\nu_{i}}\to\bar{d}$ along a subsequence. Since
$\tau_{\nu_{i}}\to 0^{+}$, we can assume that $\tau_{\nu_{i}}<t_{\min}$ for
all $i\in\mathbb{N}$, and hence Step 5 of Algorithm 1 ensures the inequality

$\displaystyle\varphi(x_{\nu_{i}}+\beta^{-1}\tau_{\nu_{i}}d_{\nu_{i}})>\varphi(x_{\nu_{i}})+\sigma\beta^{-1}\tau_{\nu_{i}}\langle
w_{\nu_{i}},d_{\nu_{i}}\rangle,\quad i\in\mathbb{N}.$ (30)

Lemma 2 gives us a constant $r>0$ such that

$\displaystyle\varphi(x_{\nu_{i}}+\beta^{-1}\tau_{\nu_{i}}d_{\nu_{i}})\leq\varphi(x_{\nu_{i}})+\beta^{-1}\tau_{\nu_{i}}\langle
w_{\nu_{i}},d_{\nu_{i}}\rangle+r\beta^{-2}\tau_{\nu_{i}}^{2}\|d_{\nu_{i}}\|^{2}$
(31)

for all $i$ sufficiently large. Combining (30), (31), and (28) tells us that

$\begin{array}[]{ll}\sigma\beta^{-1}\tau_{\nu_{i}}\langle
w_{\nu_{i}},d_{\nu_{i}}\rangle<\varphi(x_{\nu_{i}}+\beta^{-1}\tau_{\nu_{i}}d_{\nu_{i}})-\varphi(x_{\nu_{i}})\\\
\leq\beta^{-1}\tau_{\nu_{i}}\langle
w_{\nu_{i}},d_{\nu_{i}}\rangle+r\beta^{-2}\tau_{\nu_{i}}^{2}\|d_{\nu_{i}}\|^{2}\leq\beta^{-1}\tau_{\nu_{i}}\left(1-\displaystyle\frac{r}{\zeta\beta}\tau_{\nu_{i}}\right)\langle
w_{\nu_{i}},d_{\nu_{i}}\rangle\end{array}$

for large $i$. Since $\langle w_{\nu_{i}},d_{\nu_{i}}\rangle<0$ by (28), we
get that $\sigma>1-\frac{r}{\zeta\beta}\tau_{\nu_{i}}$ for such $i$, which
contradicts the choice of $\sigma\in(0,1)$ and thus verifies this claim.

Claim 2: _We have the series convergence
$\sum_{j\in\mathbb{N}}\|d_{k_{j}}\|^{2}<\infty$,
$\sum_{j\in\mathbb{N}}\|x_{k_{j}+1}-x_{k_{j}}\|^{2}<\infty$, and
$\sum_{j\in\mathbb{N}}\|w_{k_{j}}\|^{2}<\infty$._
To justify this, deduce from Step 5 of Algorithm 1 and Lemma 1 that

$\displaystyle\sum\limits_{k\in\mathbb{N}}\zeta\tau_{k}\|d_{k}\|^{2}\leq\frac{1}{\sigma}\Big{(}\varphi(x_{0})-\inf_{k\in\mathbb{N}}\varphi(x_{k})\Big{)}.$

It follows from Claim 1 that $\zeta\tau_{k_{j}}>\gamma>0$ for all
$j\in\mathbb{N}$, which yields
$\sum_{j\in\mathbb{N}}\|d_{k_{j}}\|^{2}<\infty$. On the other hand, we have
that $\|x_{k_{j}+1}-x_{k_{j}}\|=\tau_{k_{j}}\|d_{k_{j}}\|$, and again Claim 1
ensures that $\sum_{j\in\mathbb{N}}\|x_{k_{j}+1}-x_{k_{j}}\|^{2}<\infty$. To
proceed further, let $l_{2}:=\sup\\{\|d_{k_{j}}\|\;|\;\in\mathbb{N}\\}$ and
use the Lipschitz continuity of $\nabla g$ on the compact set ${\rm
cl}\\{x_{k_{j}}\;|\;{j\in\mathbb{N}}\\}\subseteq\Omega$. Employing the
subdifferential condition from (MR3823783, , Theorem 4.15) together with the
coderivative scalarization in (6), we get by the standard compactness argument
the existence of $l_{3}>0$ such that

$\displaystyle w\in\partial\langle d,\nabla
g\rangle(x_{k_{j}})=\partial^{2}g(x_{k_{j}})(d)\Longrightarrow\|w\|\leq l_{3}$

for all $j\in\mathbb{N}$ and all $d\in\mathbb{B}_{l_{2}}(0)$. Therefore, it
follows from the inclusion
$-w_{k_{j}}\in\partial^{2}g(x_{k_{j}})(d_{k_{j}})+\rho_{k_{j}}d_{k_{j}}$ that
we have

$\displaystyle\|w_{k_{j}}+\rho_{k_{j}}d_{k_{j}}\|\leq
l_{3}\|d_{k_{j}}\|\;\text{ for all large }\;j\in\mathbb{N}.$ (32)

Using finally the triangle inequality and the estimate
$\rho_{k}\leq\rho_{\max}$ leads us to the series convergence
$\sum_{j\in\mathbb{N}}\|w_{k_{j}}\|^{2}<\infty$ as stated in Claim 2.

Claim 3: _If the sequence $\\{x_{k}\\}$ is bounded, then the set of its
accumulation points is nonempty, closed and connected._
Applying Claim 2 to the sequence $\\{x_{k}\\}$, we have the _Ostrowski
condition_ $\lim_{k\to\infty}\|x_{k+1}-x_{k}\|=0$. Then, the conclusion
follows from (Ostrowski1966, , Theorem 28.1).

Claim 4: _If $x_{k_{j}}\to\bar{x}$ as $j\to\infty$, then $\bar{x}$ is a
stationary point of (1) being such that
$\varphi(\bar{x})=\inf_{k\in\mathbb{N}}\varphi(x_{k})$._
By Claim 2, we have that the sequence $w_{k_{j}}\in\partial\varphi(x_{k_{j}})$
with $w_{k_{j}}\to 0$ as $j\to\infty$. The closedness of the basic subgradient
set ensures that $0\in\partial\varphi(\bar{x})$. The second assertion of the
claim follows from the continuity of $\varphi$ at $\bar{x}\in\Omega$.

Claim 5: _If $\\{x_{k}\\}$ has an isolated accumulation point $\bar{x}$, then
the entire sequence of $x_{k}$ converges to $\bar{x}$ as $k\to\infty$, and
$\bar{x}$ is a stationary point of (1)._ Indeed, consider any subsequence
$x_{k_{j}}\to\bar{x}$. By Claim 4, $\bar{x}$ is a stationary point of (1), and
it follows from Claim 2 that $\lim_{j\to\infty}\|x_{k_{j}+1}-x_{k_{j}}\|=0$.
Then we deduce from by (MR1955649, , Proposition 8.3.10) that
$x_{k}\to\bar{x}$ as $k\to\infty$, which completes the proof of theorem.

###### Remark 7

Regarding Theorem 3.1, observe the following:

(i) If $h=0$, $g$ is of class $\mathcal{C}^{1,1}$, and $\xi>0$, then the
results of Theorem 3.1 can be found in 2021arXiv210902093D .

(ii) If $\xi\geq 0$, we can choose the regularization parameter
$\rho_{k}:=c\|w_{k}\|$ and (a varying) $\zeta:=c\|w_{k}\|$ in (23) for some
$c>0$ to verify that assertions (i) and (iii) of Theorem 3.1 still hold.
Indeed, if $\\{x_{k_{j}}\\}$ converges to some $\bar{x}$, then
$\\{w_{k_{j}}\\}$ is bounded by the Lipschitz continuity of $\varphi$. Hence
the sequence $\\{w_{k_{j}}\\}$ converges to $0$. Otherwise, there exists $M>0$
and a subsequence of $\\{w_{k_{j}}\\}$ whose norms are bounded from below by
$M$. Using the same argumentation as in the proof of Theorem 3.1 with
$\zeta=cM$, we arrive at the contradiction with $0$ being an accumulation
point of of $\\{w_{k_{j}}\\}$.

When the objective function $\varphi$ is coercive and its stationary points
are isolated, Algorithm 1 converges to a stationary point because Theorem
3.1(iii) ensures that the set of accumulation points is connected. This
property enables us to prove the convergence in some settings when even there
exist nonisolated accumulation points; see the two examples below.

###### Example 3

Consider the function $\varphi:\mathbb{R}\to\mathbb{R}$ given by

$\displaystyle\varphi(x)$
$\displaystyle:=\int_{0}^{x}t^{4}\sin\left(\frac{\pi}{t}\right)dt.$

This function is clearly $\mathcal{C}^{2}$-smooth and coercive. For any
starting point $x_{0}$, the level set
$\Omega=\\{x\;|\;\varphi(x)\leq\varphi(x_{0})\\}$ is bounded, and hence there
exists a number $\xi\in\mathbb{R}$ such that the functions $g(x):=\varphi(x)$
and $h(x):=0$ satisfy the assumptions of Theorem 3.1. Observe furthermore that
$\varphi$ is a DC function because it is $\mathcal{C}^{2}$-smooth; see, e.g.,
Oliveira_2020 ; hiriart . However, it is not possible to write its DC
decomposition with $g(x)=\varphi(x)+ax^{2}$ and $h(x)=ax^{2}$ for $a>0$, since
there exists no scalar $a>0$ such that the function $g(x)=\varphi(x)+ax^{2}$
is convex on the entire real line.

It is easy to see that the stationary points of $\varphi$ are described by
$S:=\left\\{\frac{1}{n}\;\big{|}\;n\in\mathbb{Z}\backslash\\{0\\}\right\\}\cup\\{0\\}$.
Moreover, if Algorithm 1 generates an iterative sequence $\\{x_{k}\\}$
starting from $x_{0}$, then the accumulation points form by Theorem 3.1(ii) a
nonempty, closed, and connected set $A\subseteq S$ . If $A=\\{0\\}$, the
sequence $\\{x_{k}\\}$ converges to $\bar{x}=0$. If $A$ contains any point of
the form $\bar{x}=\frac{1}{n}$, then it is an isolated point, and Theorem
3.1(iv) tells us that the entire sequence $\\{x_{k}\\}$ converges to that
point, and consequently we have $A=\\{\bar{x}\\}$.

###### Example 4

Consider the function $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ given by

$\displaystyle\varphi(x):=\sum_{i=1}^{n}\varphi_{i}(x_{i}),\;\text{ where
}\;\varphi_{i}(x_{i}):=g_{i}(x_{i})-h_{i}(x_{i})$ $\displaystyle\text{ with
}\;g_{i}(x_{i}):=\frac{1}{2}x_{i}^{2}\;\text{ and
}\;h_{i}(x_{i}):=|x_{i}|+\big{|}1-|x_{i}|\,\big{|}.$

We can easily check that the function $\varphi$ is coercive and satisfies the
assumptions of Theorem 3.1 with $g(x):=\sum_{i=1}^{n}g_{i}(x_{i})$,
$h(x):=\sum_{i=1}^{n}h_{i}(x_{i})$, and $\xi=1$. For this function, the points
in the set $\\{-2,-1,0,1,2\\}^{n}$ are critical but not stationary. Moreover,
the points in the set $\\{-2,0,2\\}^{n}$ give the global minima to the
objective function $\varphi$. Therefore, Algorithm 1 leads us to global
minimizers of $\varphi$ starting from any initial point.

$-4$$-3$$-2$$-1$$1$$2$$3$$4$$-1$$1$$2$$x$$y$ Figure 1:  Plot of the function
$\varphi_{i}$ in Example 4

The following theorem establishes convergence rates of the iterative sequences
in Algorithm 1 under some additional assumptions.

###### Theorem 3.2

Suppose in addition to the assumptions of Theorem 3.1, that $\\{x_{k}\\}$ has
an accumulation point $\bar{x}$ such that the subgradient mapping
$\partial\varphi$ is strongly metrically subregular at $(\bar{x},0)$. Then the
entire sequence $\\{x_{k}\\}$ converges to $\bar{x}$ with the Q-linear
convergence rate for $\\{\varphi(x_{k})\\}$ and the R-linear convergence rate
for $\\{x_{k}\\}$ and $\\{w_{k}\\}$. If furthermore, $\xi>0$,
$0<\zeta\leq\xi$, $\rho_{k}\to 0$, $\sigma\in(0,\frac{1}{2})$, $t_{\min}=1$,
$g$ is semismoothly differentiable at $\bar{x}$, $h$ is of class
$\mathcal{C}^{1,1}$ around $\bar{x}$, and $\mbox{\rm clm}\,\nabla
h(\bar{x})=0$, then the rate of convergence of all the sequences above is at
least Q-superlinear.

###### Proof

We split the proof of the theorem into the following two claims.

Claim 1: _The rate of convergence of $\\{\varphi(x_{k})\\}$ is at least
Q-linear, while both sequences $\\{x_{k}\\}$ and $\\{w_{k}\\}$ converge at
least R-linearly._
Observe first that it follows from the imposed strong metric subregularity of
$\partial\varphi$ that $\bar{x}$ is an isolated accumulation point, and so
$x_{k}\to\bar{x}$ as $k\to\infty$ by Theorem 3.1(iii). Further, we get from
(7) that there exists $\kappa>0$ such that

$\displaystyle\|x_{k}-\bar{x}\|\leq\kappa\|w_{k}\|\;\text{ for large
}\;k\in\mathbb{N},$ (33)

since $w_{k}\to 0$ as $k\to\infty$ by Theorem 3.1(ii). Using (32) and the
triangle inequality gives us $\ell>0$ such that $\|w_{k}\|\leq\ell\|d_{k}\|$
for sufficiently large $k\in\mathbb{N}$. Lemma 2 yields then the cost function
increment estimate

$\displaystyle\varphi(x_{k})-\varphi(\bar{x})\leq
r\|x_{k}-\bar{x}\|^{2}\;\text{ for all large }\;k\in\mathbb{N}.$ (34)

By Step 5 of Algorithm 1 and Lemma 1, we get that
$\varphi(x_{k})-\varphi(x_{k+1})\geq\sigma\zeta\tau_{k}\|d_{k}\|^{2}$ for
large $k\in\mathbb{N}$. Remembering that $\inf_{k\in\mathbb{N}}\tau_{k}>0$, we
deduce from Theorem 3.1(ii) the existence of $\eta>0$ such that

$\displaystyle\varphi(x_{k})-\varphi(\bar{x})-(\varphi(x_{k+1})-\varphi(\bar{x}))\geq\eta\|w_{k}\|^{2}$
(35)

whenever $k$ large enough. Therefore, applying (2021arXiv210902093D, , Lemma
7.2) to the sequences $\alpha_{k}:=\varphi(x_{k})-\varphi(\bar{x})$,
$\beta_{k}:=\|w_{k}\|$, and $\gamma_{k}:=\|x_{k}-\bar{x}\|$ with the positive
constants $c_{1}:=\eta$, $c_{2}:=\kappa^{-1}$, and $c_{3}:=r$, we verify the
claimed result.

Claim 2:  _Assuming that $\sigma\in(0,\frac{1}{2})$, $t_{\min}=1$, $g$ is
semismoothly differentiable at $\bar{x}$, $h$ is of class $\mathcal{C}^{1,1}$
around $\bar{x}$, and $\mbox{\rm clm}\,\nabla h(\bar{x})=0$, we have that the
rate of convergence for all the above sequences is at least Q-superlinear._
Suppose without loss of generality that $h$ is differentiable at any
$x_{k}\to\bar{x}$. It follows from the coderivative scalarization (6) and the
basic subdifferential sum rule in (MR3823783, , Theorem 2.19) valid under the
imposed assumptions that

$\displaystyle\partial^{2}g(x_{k})(d_{k})\subseteq\partial^{2}g(x_{k})(x_{k}+d_{k}-\bar{x})+\partial^{2}g(x_{k})(-x_{k}+\bar{x}).$
(36)

This yields the existence of
$z_{k}\in\partial^{2}g(x_{k})(-x_{k}+\bar{x})+\rho_{k}(-x_{k}+\bar{x})$ such
that

$\displaystyle-\nabla g(x_{k})+\nabla
h(x_{k})-z_{k}\in\partial^{2}g(x_{k})(x_{k}+d_{k}-\bar{x})+\rho_{k}(x_{k}+d_{k}-\bar{x}).$
(37)

Moreover, the $(\xi+\rho_{k})$-lower-definiteness of
$\partial^{2}g(x_{k})+\rho_{k}I$ and the Cauchy–Schwarz inequality imply that

$\displaystyle\|x_{k}+d_{k}-\bar{x}\|\leq\frac{1}{\xi+\rho_{k}}\|\nabla
g(x_{k})-\nabla h(x_{k})+z_{k}\|.$

Combining now the semismoothness of $\nabla g$ at $\bar{x}$ with the
conditions $\nabla g(\bar{x})=\nabla h(\bar{x})$ and $\mbox{\rm clm}\,\nabla
h(\bar{x})=0$ brings us to the estimates

$\displaystyle\begin{array}[]{ll}\|\nabla g(x_{k})-\nabla
h(x_{k})+z_{k}\|\leq\|\nabla g(x_{k})-\nabla
g(\bar{x})+z_{k}+\rho_{k}(x_{k}-\bar{x})\|\\\
+\rho_{k}\|x_{k}-\bar{x}\|+\|\nabla h(\bar{x})-\nabla
h(x_{k})\|=o(\|x_{k}-\bar{x}\|).\end{array}$

Then we have $\|x_{k}+d_{k}-\bar{x}\|=o(\|x_{k}-\bar{x}\|)$ and deduce
therefore from (MR1955649, , Proposition 8.3.18) and Lemma 1(i) that

$\displaystyle\varphi(x_{k}+d_{k})\leq\varphi(x_{k})+\sigma\langle\nabla\varphi(x_{k}),d_{k}\rangle.$
(38)

It follows from (38) that $x_{k+1}=x_{k}+d_{k}$ if $k$ for large $k$. Applying
(MR1955649, , Proposition 8.3.14) yields the $Q$-superlinear convergence of
$\\{x_{k}\\}$ to $\bar{x}$ as $k\to\infty$.

Finally, conditions (33)–(35) and the Lipschitz continuity of $\nabla\varphi$
around $\bar{x}$ ensure the existence of $L>0$ such that

$\displaystyle\frac{\eta}{\kappa^{2}}\|x_{k}-\bar{x}\|^{2}$
$\displaystyle\leq\varphi(x_{k})-\varphi(\bar{x})\leq r\|x_{k}-\bar{x}\|^{2},$
$\displaystyle\quad\frac{1}{\kappa}\|x_{k}-\bar{x}\|$
$\displaystyle\leq\|\nabla\varphi(x_{k})\|\leq L\|x_{k}-\bar{x}\|$

for sufficiently large $k$, and therefore we get the estimates

$\begin{array}[]{ll}\displaystyle\frac{\varphi(x_{k+1})-\varphi(\bar{x})}{\varphi(x_{k})-\varphi(\bar{x})}\leq\kappa
r\displaystyle\frac{\|x_{k+1}-\bar{x}\|^{2}}{\|x_{k}-\bar{x}\|^{2}},\\\
\quad\;\displaystyle\frac{\|\nabla\varphi(x_{k+1})\|}{\|\nabla\varphi(x_{k})\|}\leq\kappa
L\displaystyle\frac{\|x_{k+1}-\bar{x}\|}{\|x_{k}-\bar{x}\|},\end{array}$ (39)

which thus conclude the proof of the theorem.

###### Remark 8

The property of strong metric subregularity of subgradient mappings, which is
a central assumption of Theorem 3.2, has been well investigated in variational
analysis, characterized via second-order growth and coderivative type
conditions, and applied to optimization-related problems; see, e.g., ag ; dmn
; MR3823783 and the references therein.

The next theorem establishes the $Q$-superlinear and $Q$-quadratic convergence
of the sequences generated by Algorithm 1 provided that: $\xi>0$ (i.e.,
$\partial^{2}g(x)$ is $\xi$-strongly positive-definite), $\rho_{k}=0$ for all
$k\in\mathbb{N}$ (no regularization is used), $g$ is semismoothly
differentiable at the cluster point $\bar{x}$, and the function $h$ can be
expressed as the pointwise maximum of finitely many affine functions at
$\bar{x}$, i.e., when there exist
$(x^{\ast}_{i},\alpha_{i})_{i=1}^{p}\subseteq\mathbb{R}^{n}\times\mathbb{R}$
and $\varepsilon>0$ such that

$\displaystyle h(x)=\max_{i=1,\ldots,p}\left\\{\langle
x^{\ast}_{i},x\rangle+\alpha_{i}\right\\}\;\text{ for all
}\;x\in\mathbb{B}_{\varepsilon}(\bar{x}).$ (40)

###### Theorem 3.3

In addition to the assumptions of Theorem 3.1, suppose that $\xi>0$,
$0<\zeta\leq\xi$, $\sigma\in(0,\frac{1}{2})$, $t_{\min}=1$, and $\rho_{k}=0$
for all $k\in\mathbb{N}$. Suppose also that the sequence $\\{x_{k}\\}$
generated by Algorithm 1 has an accumulation point $\bar{x}$ at which $g$ is
semismoothly differentiable and $h$ can be represented in form (40). Then we
have the convergence $x_{k}\to\bar{x}$, $\varphi(x_{k})\to\varphi(\bar{x})$,
$w_{k}\to 0$, and $\nabla g(x_{k})\to\nabla g(\bar{x})$ as $k\to\infty$ with
at least $Q$-superlinear rate. If in addition $g$ is of class
$\mathcal{C}^{2,1}$ around $\bar{x}$, then the rate of convergence is at least
quadratic.

###### Proof

Observe that by (40) and (MR2191744, , Proposition 1.113) we have the
inclusion

$\displaystyle\partial(-h)(x)\subseteq\bigcup\big{\\{}-x^{\ast}_{i}\;\big{|}\;h(x)=\langle
x^{\ast}_{i},x\rangle+\alpha_{i}\big{\\}}$ (41)

for all $x$ near $\bar{x}$. The rest of the proof is split into the five
claims below.

Claim 1: _The sequence $\\{x_{k}\\}$ converges to $\bar{x}$ as $k\to\infty$._
Observe that $\bar{x}$ is an isolated accumulation point. Indeed, suppose on
the contrary that there is a sequence $\\{y_{\nu}\\}$ of accumulation points
of $\\{x_{k}\\}$ such that $y_{\nu}\to\bar{x}$ as $\nu\to\infty$ with
$y_{\nu}\neq\bar{x}$ for all $\nu\in\mathbb{N}$. Since each $y_{\nu}$ is
accumulation point of $\\{x_{k}\\}$, they are stationary points of $\varphi$.
The ${\cal C}^{1}$-smoothness of $g$ ensures that $\nabla g(y_{\nu})\to\nabla
g(\bar{x})$ as $\nu\to\infty$, and so (41) yields $\nabla
g(y_{\nu})=x_{i_{\nu}}^{\ast}$ for large $\nu\in\mathbb{N}$. Since there are
finitely many of $x_{i}^{\ast}$ in (40), we get that $\nabla g(y_{\nu})=\nabla
g(\bar{x})$ when $\nu$ is sufficiently large. Further, it follows from
(MR3823783, , Theorem 5.16) that the gradient mapping $\nabla g$ is strongly
locally maximal monotone around $\bar{x}$, i.e., there exist positive numbers
$\varepsilon$ and $r$ such that

$\displaystyle\langle\nabla g(x)-\nabla g(y),x-y\geq r\|x-y\|^{2}\;\text{ for
all }\;x,y\in\mathbb{B}_{\varepsilon}(\bar{x}).$

Putting $x:=\bar{x}$ and $y:=y_{\nu}$ in the above inequality tells us that
$\bar{x}=y_{\nu}$ for large $\nu\in\mathbb{N}$, which is a contradiction.
Applying finally Theorem 3.1(iv), we complete the proof of this claim.

Claim 2: _The sequence $\\{x_{k}\\}$ converges to $\bar{x}$ as $k\to\infty$ at
least $Q$-superlinearly._
As $x_{k}\to\bar{x}$, we have by Theorem 3.1(ii) that $w_{k}-\nabla
g(x_{k})\to-\nabla g(\bar{x})$, and so it follows from (41) that there exists
$i\in\\{1,\ldots,p\\}$ such that $h(\bar{x})=\langle
x^{\ast}_{i},\bar{x}\rangle+\alpha_{i}$, $h(x_{k})=\langle
x^{\ast}_{i},x_{k}\rangle-\alpha_{i}$ and $w_{k}-\nabla g(x_{k})=-\nabla
g(\bar{x})=-x_{i}^{\ast}$ for all $k$ sufficiently large. Define the auxiliary
function $\widehat{\varphi}:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ by

$\displaystyle\widehat{\varphi}(x):=g(x)-\langle
x^{\ast}_{i},x\rangle-\alpha_{i}$ (42)

and observe that $\widehat{\varphi}$ is $\mathcal{C}^{1,1}$ around $\bar{x}$
and semismoothly differentiable at this point. We have the equalities

$\displaystyle\varphi(x_{k})=\widehat{\varphi}(x_{k}),\;\varphi(\bar{x})=\widehat{\varphi}(\bar{x}),\;\nabla\widehat{\varphi}(x_{k})=w_{k},\;\text{
and }\;\nabla\widehat{\varphi}(\bar{x})=0$ (43)

for large $k$. It follows from
$\partial^{2}\widehat{\varphi}(x)=\partial^{2}g(x)$ that the mapping
$\partial^{2}\widehat{\varphi}(\bar{x})+\rho_{k}I$ is $(\xi+\rho_{k})$-lower-
definite. Using (36) and (37) with the replacement of $g$ by
$\widehat{\varphi}$ and taking (43) into account ensures the existence of
$z_{k}\in\partial^{2}\widehat{\varphi}(x_{k})(-x_{k}+\bar{x})+\rho_{k}(-x_{k}+\bar{x})$
satisfying the estimate

$\displaystyle\|x_{k}+d_{k}-\bar{x}\|\leq\frac{1}{\xi+\rho_{k}}\|\nabla\widehat{\varphi}(x_{k})-\nabla\widehat{\varphi}(\bar{x})+z_{k}\|.$

The triangle inequality and the semismoothness of $\nabla\widehat{\varphi}$ at
$\bar{x}$ yield

$\displaystyle\|\nabla\widehat{\varphi}(x_{k})-\nabla\widehat{\varphi}(\bar{x})+z_{k}\|$
$\displaystyle\leq\|\nabla\widehat{\varphi}(x_{k})-\nabla\widehat{\varphi}(\bar{x})+z_{k}+\rho_{k}(x_{k}-\bar{x})\|+\rho_{k}\|x_{k}-\bar{x}\|$
$\displaystyle=o(\|x_{k}-\bar{x}\|),$

which tells us that $\|x_{k}+d_{k}-\bar{x}\|=o(\|x_{k}-\bar{x}\|)$. Then it
follows from(MR1955649, , Proposition 8.3.18) and Lemma 1(i) above that

$\displaystyle\widehat{\varphi}(x_{k}+d_{k})\leq\widehat{\varphi}(x_{k})+\sigma\langle\nabla\widehat{\varphi}(x_{k}),d_{k}\rangle$
(44)

whenever $k$ is sufficiently large. Applying finally (MR1955649, , Proposition
8.3.14) verifies the claimed $Q$-superlinear convergence of $\\{x_{k}\\}$ to
$\bar{x}$.

Claim 3: _The gradient mapping of $\widehat{\varphi}$ from (42) is strongly
metrically regular around $(\bar{x},0)$ and hence strongly metrically
subregular at this point._
Using the $\xi$-lower-definiteness of $\partial^{2}\widehat{\varphi}(\bar{x})$
and the pointbased coderivative characterization of strong local maximal
monotonicity given in (MR3823783, , Theorem 5.16), we verify this property for
$\nabla\widehat{\varphi}$ around $\bar{x}$. Then (MR3823783, , Corollary 5.15)
ensures that $\nabla\widehat{\varphi}$ is strongly metrically regular around
$(\bar{x},0)$.

Claim 4: _The sequences $\\{\varphi(x_{k})\\}$, $\\{w_{k}\\}$, and $\\{\nabla
g(x_{k})\\}$ converge at least Q-superlinearly to $\varphi(\bar{x})$, $0$, and
$\nabla g(\bar{x})$, respectively._
It follows from the estimates in (39), with the replacement of $\varphi$ by
$\widehat{\varphi}$ and with taking into account that
$\widehat{\varphi}(x_{k})-\widehat{\varphi}(\bar{x})=\varphi(x_{k})-\varphi(\bar{x})$
and $\nabla\widehat{\varphi}(x_{k})=w_{k}$ due to (43), that there exist
constants $\alpha_{1},\alpha_{2}>0$ such that

$\displaystyle\frac{\varphi(x_{k+1})-\varphi(\bar{x})}{\varphi(x_{k})-\varphi(\bar{x})}$
$\displaystyle\leq\alpha_{1}\frac{\|x_{k+1}-\bar{x}\|^{2}}{\|x_{k}-\bar{x}\|^{2}}$
$\displaystyle\frac{\|w_{k+1}\|}{\|w_{k}\|}$
$\displaystyle\leq\alpha_{2}\frac{\|x_{k+1}-\bar{x}\|}{\|x_{k}-\bar{x}\|}$

provided that $k$ is sufficiently large. Recalling that $w_{k}-\nabla
g(x_{k})=-\nabla g(\bar{x})$ for large $k$ completes the proof of the claim.

Claim 5: _If $g$ is of class $\mathcal{C}^{2,1}$ around $\bar{x}$, then the
rate of convergence of the sequences above is at least quadratic._
It is easy to see that the assumed $\mathcal{C}^{2,1}$ property of $g$ yields
this property of $\widehat{\varphi}$ around $\bar{x}$. Using estimate (44), we
deduce this claim from the quadratic convergence of the classical Newton
method; see, e.g., (Aragon2019, , Theorem 5.18) and (MR3289054, , Theorem
2.15). This therefore completes the proof of the theorem.

###### Remark 9

Concerning Theorem 3.3, observe the following:

(i) It is important to emphasize that the performance of Algorithm 1 revealed
in Theorem 3.3 is mainly due to the usage of the basic subdifferential of the
function $-h$ in contrast to that of $h$, which is calculated as

$\partial
h(x)=\text{co}\left(\bigcup\left\\{x^{\ast}_{i}\;\bigg{|}\;h(x)=\langle
x^{\ast}_{i},x\rangle+\alpha_{i}\right\\}\right)$ (45)

by (MR2191744, , Theorem 3.46). We can see from the proof of Theorem 3.3 that
it fails if the evaluation of $\partial(-h)(x)$ in (41) is replaced by the one
of $\partial h(x)$ in (45).

(ii) The main assumptions of Theorem 3.3 do not imply the smoothness of
$\varphi$ at stationary points. For instance, consider the nonconvex function
$\varphi:\mathbb{R}^{n}\to\mathbb{R}$ defined as in Example 4 but letting now
$h_{i}(x_{i}):=|x_{i}|+|1-x_{i}|$. The function $\varphi$ satisfies the
assumptions of Theorem 3.3 at any of its stationary points $\\{-2,0,2\\}^{n}$,
but $\varphi$ is not differentiable at $\bar{x}=0$; see Figure 2.

$-5$$-4$$-3$$-2$$-1$$1$$2$$3$$4$$5$$-3$$-2$$-1$$1$$2$$x$$y$ Figure 2:  Plot of
function $\varphi_{i}(x)=\frac{1}{2}x^{2}-|x|-|1-x|$ in Remark 9

(iii) The functions $\varphi$, $g$, and $h$ in Example 4 satisfy the
assumptions of Theorem 3.3. Therefore, the convergence of the sequences
generated by Algorithm 1 is at least quadratic.

##  4 Convergence Rates under the Kurdyka–Łojasiewicz Property

In this section, we verify the global convergence of Algorithm 1 and establish
convergence rates in the general setting of Theorem 3.1 without additional
assumptions of Theorems 3.2 and 3.3 while supposing instead that the cost
function $\varphi$ satisfies the Kurdyka–Łojasiewicz property. Recall that the
_Kurdyka–Łojasiewicz property_ holds for $\varphi$ at $\bar{x}$ if there exist
$\eta>0$ and a continuous concave function $\psi:[0,\eta]\to[0,\infty)$ with
$\psi(0)=0$ such that $\psi$ is $\mathcal{C}^{1}$-smooth on $(0,\eta)$ with
the strictly positive derivative $\psi^{\prime}$ and that

$\displaystyle\psi^{\prime}\big{(}\varphi(x)-\varphi(\bar{x})\big{)}\,{\rm
dist}\big{(}0;\partial\varphi(x)\big{)}\geq 1$ (46)

for all $x\in\mathbb{B}_{\eta}(\bar{x})$ with
$\varphi(\bar{x})<\varphi(x)<\varphi(\bar{x})+\eta$, where ${\rm
dist}(\cdot;\Omega)$ stands for the distance function of a set $\Omega$.

The first theorem of this section establishes the global convergence of
iterative sequence generated by Algorithm 1 to a stationary point of (1).

###### Theorem 4.1

In addition to the assumptions of Theorem 3.1, suppose that the iterative
sequence $\\{x_{k}\\}$ generated by Algorithm 1 has an accumulation point
$\bar{x}$ at which the Kurdyka–Łojasiewicz property (46) is satisfied. Then
$\\{x_{k}\\}$ converges $\bar{x}$ as $k\to\infty$, which is a stationary point
of problem (1).

###### Proof

If Algorithm 1 stops after a finite number of iterations, there is nothing to
prove. Due to the decreasing property of $\\{\varphi(x_{k})\\}$ from Theorem
3.1(i), we can assume that $\varphi(x_{k})>\varphi(x_{k+1})$ for all
$k\in\mathbb{N}$. Let $\bar{x}$ be the accumulation point of $\\{x_{k}\\}$
where $\varphi$ satisfies the Kurdyka–Łojasiewicz inequality (46), which by
Theorem 3.1 is a stationary point of problem (1). Since $\varphi$ is
continuous, we have that
$\varphi(\bar{x})=\inf_{k\in\mathbb{N}}\varphi(x_{k})$. Taking the constant
$\eta>0$ and the function $\psi$ from (46) and remembering that $g$ is of
class $\mathcal{C}^{1,1}$ around $\bar{x}$, suppose without loss of generality
that $\nabla g$ is Lipschitz continuous on $\mathbb{B}_{2\eta}(\bar{x})$ with
modulus $\kappa$. Let $k_{0}\in\mathbb{N}$ be such that
$x_{k_{0}}\in\mathbb{B}_{\eta/2}(\bar{x})$ and that

$\displaystyle\varphi(\bar{x})<\varphi(x_{k})<\varphi(\bar{x})+\eta,\quad\frac{\kappa+\rho_{\max}}{\sigma\zeta}\psi\big{(}\varphi(x_{k}-\varphi(\bar{x})\big{)}<\eta/2$
(47)

for all $k\geq k_{0}$, where $\sigma\in(0,1)$, $\zeta>0$, and $\rho_{\max}>0$
are the constants of Algorithm 1. The rest of the proof is split into the
following three steps.

Claim 1: _Let $k\geq k_{0}$ be such that $x_{k}\in\mathbb{B}_{\eta}(\bar{x})$.
Then we have the estimate_

$\displaystyle\|x_{k}-x_{k+1}\|\leq\frac{\kappa+\rho_{k}}{\sigma\zeta}\big{(}\psi(\varphi(x_{k})-\varphi(\bar{x})\big{)}-\psi\big{(}\varphi(x_{k+1})-\varphi(\bar{x})\big{)}\big{)}.$
(48)

Indeed, it follows from (6), (23), and (MR3823783, , Theorem 1.22) that

$\begin{array}[]{ll}{\rm
dist}(0;\partial\varphi(x_{k})\big{)}&\leq\|w_{k}\|\leq\|w_{k}+\rho_{k}d_{k}\|+\rho_{k}\|d_{k}\|\\\
&\leq(\kappa+\rho_{k})\|d_{k}\|=\displaystyle\frac{\kappa+\rho_{k}}{\tau_{k}}\|x_{k+1}-x_{k}\|.\end{array}$
(49)

Then using Step 5 of Algorithm 1, Lemma 1, the Kurdyka–Łojasiewicz inequality
(46), the concavity of $\psi$, and estimate (49) gives us

$\displaystyle\|x_{k}-$ $\displaystyle
x_{k+1}\|^{2}=\tau^{2}_{k}\|d_{k}\|^{2}\leq\frac{\tau_{k}}{\sigma\zeta}\big{(}\varphi(x_{k})-\varphi(x_{k+1})\big{)}$
$\displaystyle\leq\frac{\tau_{k}}{\sigma\zeta}{\rm
dist}\big{(}0;\partial\varphi(x_{k})\big{)}\,\psi^{\prime}\big{(}\varphi(x_{k})-\varphi(\bar{x})\big{)}\big{(}\varphi(x_{k})-\varphi(x_{k+1})\big{)}$
$\displaystyle\leq\frac{\tau_{k}}{\sigma\zeta}{\rm
dist}\big{(}0;\partial\varphi(x_{k})\big{)}\big{(}\psi(\varphi(x_{k})-\varphi(\bar{x})\big{)}-\psi(\varphi(x_{k+1})-\varphi(\bar{x})\big{)}\big{)}$
$\displaystyle\leq\frac{\kappa+\rho_{k}}{\sigma\zeta}\|x_{k+1}-x_{k}\|\big{(}\psi\big{(}\varphi(x_{k})-\varphi(\bar{x})\big{)}-\psi\big{(}\varphi(x_{k+1})-\varphi(\bar{x})\big{)}\big{)},$

which therefore verifies the claimed inequality (48).

Claim 2: _For every $k\geq k_{0}$, we have the inclusion
$x_{k}\in\mathbb{B}_{\eta}(\bar{x})$._
Suppose on the contrary that there exists $k>k_{0}$ with
$x_{k}\notin\mathbb{B}_{\eta}(\bar{x})$ and define
$\bar{k}:=\min\left\\{k>k_{o}\;\big{|}\;x_{k}\notin\mathbb{B}_{\eta}(\bar{x})\right\\}$.
Since for $k\in\\{k_{0},\ldots,\bar{k}-1\\}$ the estimate in (48) is
satisfied, we get by using (47) that

$\displaystyle\|x_{\bar{k}}-\bar{x}\|$
$\displaystyle\leq\|x_{k_{0}}-\bar{x}\|+\sum_{k=k_{0}}^{\bar{k}-1}\|x_{k}-x_{k+1}\|$
$\displaystyle\leq\|x_{k_{0}}-\bar{x}\|+\frac{\kappa+\rho_{\max}}{\sigma\zeta}\sum_{k=k_{0}}^{\bar{k}-1}\big{(}\psi\big{(}\varphi(x_{k})-\varphi(\bar{x})\big{)}-\psi\big{(}\varphi(x_{k+1})-\varphi(\bar{x})\big{)}\big{)}$
$\displaystyle\leq\|x_{k_{0}}-\bar{x}\|+\frac{\kappa+\rho_{\max}}{\sigma\zeta}\psi\big{(}\varphi(x_{k_{0}})-\varphi(\bar{x})\big{)}\leq\eta,$

which contradicts our assumption and thus verifies this claim.

Claim 3: _We have that $\sum_{k=1}^{\infty}\|x_{k}-x_{k+1}\|<\infty$, and
consequently the sequence $\\{x_{k}\\}$ converges to $\bar{x}$ as
$k\to\infty$._
It follows from Claim 1 and Claim 2 that (48) holds for all $k\geq k_{0}$.
Thus

$\displaystyle\sum_{k=1}^{\infty}\|x_{k}-x_{k+1}\|$
$\displaystyle\leq\sum_{k=1}^{k_{0}-1}\|x_{k}-x_{k+1}\|+\sum_{k=k_{0}}^{\infty}\|x_{k}-x_{k+1}\|$
$\displaystyle\leq\sum_{k=1}^{k_{0}-1}\|x_{k}-x_{k+1}\|+\frac{\kappa+\rho_{\max}}{\sigma\zeta}\psi\big{(}\varphi(x_{k_{0}})-\varphi(\bar{x})\big{)}<\infty,$

which therefore completes the proof of the theorem.

The next theorem establishes convergence rates for iterative sequence
$\\{x_{k}\\}$ in Algorithm 1 provided that the function $\psi$ in (46) is
selected in a special way. Since the proof while using Theorem 4.1, is similar
to the corresponding one from (MR4078808, , Theorem 4.9) in a different
setting, it is omitted.

###### Theorem 4.2

In addition to the assumptions of Theorem 4.1, suppose that the
Kurdyka–Łojasiewicz property (46) holds at the accumulation point $\bar{x}$
with $\psi(t):=Mt^{1-\theta}$ for some $M>0$ and $\theta\in[0,1)$. The
following assertions hold:

(i) If $\theta=0$, then the sequence $\\{x_{k}\\}$ converges in a finite
number of steps.

(ii) If $\theta\in(0,1/2]$, then the sequence $\\{x_{k}\\}$ converges at least
linearly.

(iii) If $\theta\in(1/2,1)$, then there exist $\mu>0$ and $k_{0}\in\mathbb{N}$
such that

$\displaystyle\|x_{k}-\bar{x}\|\leq\mu k^{-\frac{1-\theta}{2\theta-1}}\;\text{
for all }\;k\geq k_{0}.$

###### Remark 10

Together with our main Algorithm 1, we can consider its modification with the
replacement of $\partial(-h)(x_{k})$ by $-\partial h(x_{k})$. In this case,
the most appropriate version of the Kurdyka–Łojasiewicz inequality (46),
ensuring the fulfillment the corresponding versions of Theorem 4.1 and 4.2, is
the one

$\displaystyle\psi\big{(}\varphi(x)-\varphi(\bar{x})\big{)}\,{\rm
dist}\big{(}0;\partial^{0}\varphi(x)\big{)}\geq 1$

expressed in terms of the symmetric subdifferential $\partial^{0}\varphi(x)$
from (14). Note that the latter is surely satisfied where the symmetric
subdifferential is replaced by the generalized gradient
$\overline{\partial}\varphi(x)$, which is the convex hull of
$\partial^{0}\varphi(x)$.

##  5 Applications to Structured Constrained Optimization

In this section, we present implementations and specifications of our main
RCSN Algorithm 1 for two structured classes of optimization problems. The
first class contains functions represented as sums of two nonconvex functions
one of which is smooth, while the other is extended-real-valued. The second
class concerns minimization of smooth functions over closed constraint sets.

###  5.1 Minimization of Structured Sums

Here we consider the following class of structured optimization problems:

$\min_{x\in\mathbb{R}^{n}}\varphi(x):=f(x)+\psi(x),$ (50)

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is of class $\mathcal{C}^{2,1}$ with the
$L_{f}$-Lipschitzian gradient, and where
$\psi:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is an extended-real-valued prox-
bounded function with the threshold $\lambda_{\psi}>0$. When both functions
$f$ and $\psi$ are convex, problems of type (50) have been largely studied
under the name of “convex composite optimization” emphasizing the fact that
$f$ and $\psi$ are of completely different structures. In our case, we do not
impose any convexity of $f,\psi$ and prefer to label (50) as minimization of
structured sums to avoid any confusions with optimization of function
compositions, which are typically used in major models of variational analysis
and constrained optimization; see, e.g., MR1491362 .

In contrast to the original class of unconstrained problems of difference
programming (1), the structured sum optimization form (50) covers optimization
problems with constraints given by $x\in\mbox{\rm dom}\,\psi$. Nevertheless,
we show in what follows that the general class of problem (50) can be reduced
under the assumptions imposed above to the difference form (1) satisfying the
conditions for the required performance of Algorithm 1.

This is done by using an extended notion of envelopes introduced by Patrinos
and Bemporad in Patrinos2013 , which is now commonly referred as the _forward-
backward envelope_ ; see, e.g., MR3845278 .

###### Definition 4

Given $\varphi=f+\psi$ and $\lambda>0$, the _forward-backward envelope_ (FBE)
of the function $\varphi$ with the parameter $\lambda$ is defined by

$\displaystyle\varphi_{\lambda}(x):=\inf_{z\in\mathbb{R}^{n}}\Big{\\{}f(x)+\langle\nabla
f(x),z-x\rangle+\psi(z)+\frac{1}{2\lambda}\|z-x\|^{2}\Big{\\}}.$ (51)

Remembering the constructions of the Moreau envelope (15) and the Asplund
function (17) allows us to represent $\varphi_{\lambda}$ for every
$\lambda\in(0,\lambda_{\psi})$ as:

$\begin{array}[]{ll}\varphi_{\lambda}(x)&=f(x)-\displaystyle\frac{\lambda}{2}\|\nabla
f(x)\|^{2}+{\mathtt{e}}_{\lambda}\psi\big{(}x-\lambda\nabla f(x)\big{)}\\\
&=\displaystyle f(x)+\frac{1}{2\lambda}\|x\|^{2}-\langle\nabla
f(x),x\rangle-{{\mathtt{A}}_{\lambda}\psi}\big{(}x-\lambda\nabla
f(x)\big{)}.\end{array}$ (52)

###### Remark 11

It is not difficult to show that whenever $\nabla f$ is $L_{f}$-Lipschitz on
$\mathbb{R}^{n}$ and $\lambda\in(0,\frac{1}{L_{f}})$, the optimal values in
problems (1) and (51) are the same

$\displaystyle\inf_{x\in\mathbb{R}^{n}}\varphi_{\lambda}(x)=\inf_{x\in\mathbb{R}^{n}}\varphi(x).$
(53)

Indeed, the inequality “$\leq$” in (53) follows directly from the definition
of $\varphi_{\lambda}$. The reverse inequality in (53) is obtained by

$\begin{array}[]{ll}\displaystyle\inf_{x\in\mathbb{R}^{n}}\varphi_{\lambda}(x)=\displaystyle\inf_{x\in\mathbb{R}^{n}}\displaystyle\inf_{z\in\mathbb{R}^{n}}\Big{\\{}f(x)+\langle\nabla
f(x),z-x\rangle+\psi(z)+\displaystyle\frac{1}{2\lambda}\|z-x\|^{2}\Big{\\}}\\\
\geq\displaystyle\inf_{x\in\mathbb{R}^{n}}\displaystyle\inf_{z\in\mathbb{R}^{n}}\Big{\\{}f(z)-\frac{L_{f}}{2}\|z-x\|^{2}+\psi(z)+\displaystyle\frac{1}{2\lambda}\|z-x\|^{2}\Big{\\}}\\\
=\displaystyle\inf_{z\in\mathbb{R}^{n}}\displaystyle\inf_{x\in\mathbb{R}^{n}}\Big{\\{}f(z)+\psi(z)+\Big{(}\frac{1}{2\lambda}-\displaystyle\frac{L_{f}}{2}\Big{)}\|z-x\|^{2}\Big{\\}}=\displaystyle\inf_{z\in\mathbb{R}^{n}}\varphi(x).\end{array}$

Moreover, (53) does not hold if $\nabla f$ is not Lipschitz continuous on
$\mathbb{R}^{n}$. Indeed, consider $f(x):=\frac{1}{4}x^{4}$ and $\psi:=0$.
Then we have $\inf_{x\in\mathbb{R}^{n}}\varphi(x)=0$ while
$\varphi_{\lambda}(x)=\frac{1}{4}x^{4}-\frac{\lambda}{2}x^{6}$, which yields
$\inf_{x\in\mathbb{R}^{n}}\varphi_{\lambda}(x)=-\infty$, and so (53) fails.

The next theorem shows that FBE (51) can be written as the difference of a
$\mathcal{C}^{1,1}$ function and a Lipschitzian prox-regular function.
Furthermore, it establishes relationships between minimizers and critical
points of $\varphi$ and $\varphi_{\lambda}$.

###### Theorem 5.1

Let $\varphi=f+\psi$, where $f$ is of class $\mathcal{C}^{2,1}$ and where
$\psi$ is prox-bounded with threshold $\lambda_{\psi}>0$. Then for any
$\lambda\in(0,\lambda_{\psi})$, we have the inclusion

$\partial\varphi_{\lambda}(x)\subseteq\lambda^{-1}\big{(}I-\lambda\nabla^{2}f(x)\big{)}\big{(}x-{\mathtt{Prox}}_{\lambda\psi}\big{(}x-\lambda\nabla
f(x)\big{)}\big{)}.$ (54)

Furthermore, the following assertions are satisfied:

(i) If $x\in\mathbb{R}^{n}$ is a stationary point of $\varphi_{\lambda}$, then
$0\in\widehat{\partial}\varphi(x)$ provided that the matrix
$I-\lambda\nabla^{2}f(x)$ is nonsingular.

(ii) The FBE (51) can be written as $\varphi_{\lambda}=g-h$, where
$g(x):=f(x)+\frac{1}{2\lambda}\|x\|^{2}$ is of class $\mathcal{C}^{2,1}$, and
where $h(x):=\langle\nabla
f(x),x\rangle+{{\mathtt{A}}_{\lambda}\psi}(x-\lambda\nabla f(x))$ is locally
Lipschitzian and prox-regular on $\mathbb{R}^{n}$. Moreover, $\nabla^{2}g(x)$
is $\xi$-lower-definite for all $x\in\mathbb{R}^{n}$ with
$\xi:=\frac{1}{\lambda}-L_{f}$.

(iii) If $\psi:=\delta_{C}$ for a closed set $C$, then
$\partial(-{{\mathtt{A}}_{\lambda}\psi})=-\frac{1}{\lambda}\mathtt{P}_{C}$,
where $\mathtt{P}_{C}$ denotes the $($generally set-valued$)$ projection
operator onto $C$. In this case, inclusion (54) holds as an equality.

(iv) If both $f$ and $\psi$ are convex, we have that $\varphi_{\lambda}=g-h$,
where $g(x):=f(x)+{\mathtt{e}}_{\lambda}\psi(x-\lambda\nabla f(x))$ and
$h(x):=\frac{\lambda}{2}\|\nabla f(x)\|^{2}$ are of class $\mathcal{C}^{1,1}$
$($and hence prox-regular$)$ on $\mathbb{R}^{n}$, and that

$\displaystyle\big{\\{}x\in\mathbb{R}^{n}\;\big{|}\;\nabla\varphi_{\lambda}(x)=0\big{\\}}=\big{\\{}x\in\mathbb{R}^{n}\;\big{|}\;0\in\partial\varphi(x)\big{\\}}$
(55)

provided that $I-\lambda\nabla^{2}f(x)$ is nonsingular at any stationary point
of $\varphi_{\lambda}$.

###### Proof

Observe that inclusion (54) follows directly by applying the basic
subdifferential sum and chain rules from (MR3823783, , Theorem 2.19 and
Corollary 4.6), respectively, the first representation of $\varphi_{\lambda}$
in (52) with taking into account the results of Lemma 2. Now we pick any
stationary point $x\in\mathbb{R}^{n}$ of the FBE $\varphi_{\lambda}$ and then
deduce from $0\in\partial\varphi_{\lambda}(x)$ and (54) that

$x\in{\mathtt{Prox}}_{\lambda\psi}\big{(}x-\lambda\nabla f(x))\big{)},$

which readily implies that $0\in\nabla
f(x)+\widehat{\partial}\psi(x)=\widehat{\partial}\varphi(x)$ and thus verifies
(i). Assertion (ii) follows directly from Proposition 2 and the smoothness of
$f$.

To prove (iii), we need to verify the reverse inclusion “$\supseteq$” in (19),
for which it suffices to show that the inclusion $v\in\mathtt{P}_{C}(x)$
yields $v\not\in{\rm co}(\mathtt{P}_{C}(x)\setminus\\{v\\})$. On the contrary,
if $v\in\mathtt{P}_{C}(x)\cap{\rm co}(\mathtt{P}_{C}(x)\setminus\\{v\\})$,
then there exist $c_{1},\ldots,c_{m}\in P_{C}(x)\setminus\\{v\\}$ and
$\mu_{1},\ldots,\mu_{m}\in(0,1)$ such that $v=\sum_{i=1}^{m}\mu_{i}c_{i}$ with
$\sum_{i=1}^{m}\mu_{i}=1$. By definition of the projection, we get the
equalities

$|c_{1}-x\|^{2}=\ldots=\|c_{m}-x\|^{2}=\|v-x\|^{2}=\Big{\|}\sum_{i=1}^{m}\mu_{i}(c_{i}-x)\Big{\|}^{2},$

which contradict the strict convexity of $\|\cdot\|^{2}$ and thus verifies
(iii).

The first statement in (iv) follows from the differentiability of $f$ and of
the Moreau envelope ${\mathtt{e}}_{\lambda}\psi$ by (MR1491362, , Theorem
2.26). Further, the inclusion “$\subseteq$” in (55) is a consequence of (i).
To justify the reverse inclusion in (55), observe that any $x$ satisfying
$0\in\partial\varphi(x)$ is a global minimizer of the convex function
$\varphi$, and so $x={\mathtt{Prox}}_{\lambda\psi}(x-\lambda\nabla f(x))$. The
differentiability of $\varphi_{\lambda}$ and (54) (which holds as an equality
in this case) tells us that $\nabla\varphi_{\lambda}(x)=0$, and thus (55)
holds. This completes the proof of the theorem.

###### Remark 12

Based on Theorem 5.1(ii), it is not hard to show that the FBE function
$\varphi_{\lambda}$ can be represented as a difference of convex functions.
Indeed, since ${{\mathtt{A}}_{\lambda}\psi}$ is a locally Lipschitzian and
prox-regular function, we have by (MR2101873, , Corollary 3.12) that $h$ is a
lower-$\mathcal{C}^{2}$ function, and hence by (MR1491362, , Theorem 10.33),
it is locally a DC function. Similarly, $g$ being a $\mathcal{C}^{2}$ function
is a DC function, so the difference $\varphi=g-h$ is also a DC function.
However, it is difficult to determine for numerical purposes what is an
appropriate representation of $\varphi$ as a difference of convex functions.
Moreover, such a representation of the objective in terms of convex functions
may generate some theoretical and algorithmic challenges as demonstrated below
in Example 5.

###  5.2 Nonconvex Optimization with Geometric Constraints

This subsection addresses the following problem of constrained optimization
with explicit geometric constraints given by:

$\mbox{minimize }\;f(x)\;\mbox{ subject to }\;x\in C,$ (56)

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is of class $\mathcal{C}^{2,1}$, and
where $C\subseteq\mathbb{R}^{n}$ is an arbitrary closed set. Due to the lack
of convexity, most of the available algorithms in the literature are not able
to directly handle this problem. Nevertheless, Theorem 5.1 provides an
effective machinery allowing us to reduce (56) to an optimization problem that
can be solved by using our developments. Indeed, define
$\psi(x):=\delta_{C}(x)$ and observe that $\psi$ is prox-regular with
threshold $\lambda_{\psi}=\infty$. In this setting, FBE (51) reduces to the
formula

$\varphi_{\lambda}(x)=f(x)-\frac{\lambda}{2}\|\nabla
f(x)\|^{2}+\frac{1}{2\lambda}{\rm dist}^{2}\big{(}x-\lambda\nabla
f(x);C\big{)}.$

Furthermore, it follows from Theorem 5.1(iii) that

$\displaystyle\partial\varphi_{\lambda}(x)=\lambda^{-1}\big{(}I-\lambda\nabla^{2}f(x)\big{)}\big{(}x-\mathtt{P}_{C}\big{(}x-\lambda\nabla
f(x)\big{)}\big{)}.$

Based on Theorem 5.1, we deduce from Algorithm 1 with $\rho_{k}=0$ its
following version to solve the constrained problem (56).

1:$x_{0}\in\mathbb{R}^{n}$, $\beta\in(0,1)$, $t_{\min}>0$ and
$\sigma\in(0,1)$.

2:for $k=0,1,\ldots$ do

3: Take
$w_{k}\in\left(\lambda^{-1}I-\nabla^{2}f(x_{k})\right)\big{(}x_{k}-\mathtt{P}_{C}\big{(}x-\lambda\nabla
f(x_{k})\big{)}\big{)}$.

4: If $w_{k}=0$, STOP and return $x_{k}$. Otherwise set $d_{k}$ as the
solution to the linear system $(\nabla^{2}f(x_{k})+\lambda^{-1}I)d_{k}=w_{k}$.

5: Choose any $\overline{\tau}_{k}\geq t_{\min}$. Set
$\overline{\tau}_{k}:=\tau_{k}$.

6: while
$\varphi_{\lambda}(x_{k}+\tau_{k}d_{k})>\varphi_{\lambda}(x_{k})+\sigma\tau_{k}\langle\nabla
w_{k},d_{k}\rangle$ do

7: $\tau_{k}=\beta\tau_{k}$.

8: end while

9: Set $x_{k+1}:=x_{k}+\tau_{k}d_{k}$.

10:end for

Algorithm 2 Projected-like Newton algorithm for constrained optimization

To the best of our knowledge, Algorithm 2 is new even for the case of convex
constraint sets $C$. All the results obtained for Algorithm 1 in Sections 3
and 4 can be specified for Algorithm 2 to solve problem (56). For brevity, we
present just the following direct consequence of Theorem 3.1.

###### Corollary 1

Considering problem (56), suppose that $f:\mathbb{R}^{n}\to\mathbb{R}$ is of
class $\mathcal{C}^{2,1}$, that $C\subset\mathbb{R}^{n}$ is closed, and that
$\inf_{x\in C}f(x)>-\infty$. Pick an initial point $x_{0}\in\mathbb{R}^{n}$
and a parameter $\lambda\in(0,\frac{1}{L_{f}})$. Then Algorithm 2 either stops
at a point $x$ such that $0\in\nabla f(x)+\widehat{N}_{C}(x)$, or generates
infinite sequences $\\{x_{k}\\}$, $\\{\varphi_{\lambda}(x_{k})\\}$,
$\\{w_{k}\\}$, $\\{d_{k}\\}$, and $\\{\tau_{k}\\}$ satisfying the assertions:

(i) The sequence $\\{\varphi_{\lambda}(x_{k})\\}$ monotonically decreases and
converges.

(ii) If $\\{x_{k_{j}}\\}$ is a bounded subsequence of $\\{x_{k}\\}$, then
$\inf_{j\in\mathbb{N}}\tau_{k_{j}}>0$ and

$\displaystyle\sum\limits_{j\in\mathbb{N}}\|d_{k_{j}}\|^{2}<\infty,\;\sum\limits_{j\in\mathbb{N}}\|x_{k_{j}+1}-x_{k_{j}}\|^{2}<\infty,\;\sum\limits_{j\in\mathbb{N}}\|w_{k_{j}}\|^{2}<\infty.$

If, in particular, the entire sequence $\\{x_{k}\\}$ is bounded, then the set
of its accumulation points is nonempty, closed, and connected.

(iii) If $x_{k_{j}}\to\bar{x}$ as $j\to\infty$, then $0\in\nabla
f(\bar{x})+\widehat{N}_{C}(\bar{x})$ and the equality
$\varphi_{\lambda}(\bar{x})=\inf_{k\in\mathbb{N}}\varphi_{\lambda}(x_{k})$
holds.

(iv) If the sequence $\\{x_{k}\\}$ has an isolated accumulation point
$\bar{x}$, then it converges to $\bar{x}$ as $k\to\infty$, and we have
$0\in\nabla f(\bar{x})+\widehat{N}_{C}(\bar{x})$.

The next example illustrates our approach to solve (56) via Algorithm 2 in
contrast to algorithms of the DC type.

###### Example 5

Consider the minimization of a quadratic function over a closed (possibly
nonconvex) set $C$:

$\displaystyle\mbox{minimize }\;\frac{1}{2}x^{T}Qx+b^{T}x\;\text{ subject to
}\;x\in C,$ (57)

where $Q$ is a symmetric matrix, and where $b\in\mathbb{R}^{n}$. In this
setting, FBE (51) can be written as $\varphi_{\lambda}(x)=g(x)-h(x)$ with

$\begin{array}[]{ll}g(x)&:=\displaystyle\frac{1}{2}x^{T}\big{(}Q+\lambda^{-1}I\big{)}x+b^{T}x,\\\
h(x)&:=x^{T}Qx+b^{T}x+{{\mathtt{A}}_{\lambda}\psi}\big{(}(I-\lambda
Q)x-\lambda b\big{)}.\end{array}$ (58)

Our method does not require a DC decomposition of the objective function
$\varphi_{\lambda}$. Indeed the function $h$ in (58) is generally nonconvex.
Specifically, consider $Q=\begin{bmatrix}0&-1\\\ -1&0\end{bmatrix}$,
$b=(0,0)^{T}$, and $C$ being the unit sphere centered at the origin. Then $g$
in (58) is strongly convex for any $\lambda\in(0,1)$, while $h$ therein is not
convex whenever $\lambda>0$. More precisely, in this case we have

$h(x_{1},x_{2})=-2x_{1}x_{2}+{{\mathtt{A}}_{\lambda}\psi}(x_{1}+\lambda
x_{2},\lambda x_{1}+x_{2})\;\mbox{ with}$
$\displaystyle{{\mathtt{A}}_{\lambda}\psi}(x)=\frac{1}{2\lambda}\left(\|x\|^{2}-d_{C}^{2}(x)\right)=\frac{1}{2\lambda}\left(\|x\|^{2}-(\|x\|-1)^{2}\right)=\frac{1}{2\lambda}\left(2\|x\|-1\right).$

This tells us, in particular, that

$h(-1/2,-1/2)-\frac{1}{2}h(-1,-1)-\frac{1}{2}h(0,0)=\frac{1}{2},$

and thus $h$ is not convex regardless of the value of $\lambda$; see Figure 3.

Figure 3: Contour plot of the functions $f$, $\varphi_{\lambda}$, $g$ and $h$
in (58) with $\lambda=0.9$

##  6 Further Applications and Numerical Experiments

In this section, we demonstrate the performance of Algorithm 1 and Algorithm 2
in two different problems. The first problem is smooth and arises from the
study of system biochemical reactions. It can be successfully tackled with
DCA-like algorithms, but they require to solve subproblems whose solutions
cannot be analytically computed and are thus time-consuming. This is in
contrast to Algorithm 1, which only requires solving the linear equation (23)
at each iteration. The second problem is nonsmooth and consists of minimizing
a quadratic function under both convex and nonconvex constraints. Employing
FBE (51) and Theorem 5.1, these two problems can be attacked by using DCA,
BDCA, and Algorithm 2.

Both Algorithms 1 and 2 have complete freedom in the choice of the initial
value of the stepsizes $\overline{\tau}_{k}$ in Step 4, as long as they are
bounded from below by a positive constant $t_{\min}$, while the choice of
$\overline{\tau}_{k}$ totally determines the performance of the algorithms. On
the one hand, a small value would permit the stepsize to get easily accepted
in Step 5, but it would imply little progress in the iteration and (likely) in
the reduction of the objective function, probably making it more prone to
stagnate at local minima. On the other hand, we would expect a large value to
ameliorate these issues, while it could result in a significant waste of time
in the linesearch Steps 5-7 of both algorithms.

Therefore, it makes sense to consider a choice which sets the trial stepsize
$\overline{\tau}_{k}$ depending on the stepsize $\tau_{k-1}$ accepted in the
previous iteration, perhaps increasing it if no reduction of the stepsize was
needed. This technique was introduced in (MR4078808, , Section 5) under the
name of _Self-adaptive trial stepsize_ , and it was shown there that this
accelerates the performance of BDCA in practice. A similar idea is behind the
so-called _two-way backtracking_ linesearch, which was recently proposed in
Truong2021 for the gradient descent method, showing good numerical results on
deep neural networks. In contrast to BDCA, our theoretical results require
$t_{\min}$ to be strictly positive, so the technique should be slightly
adapted as shown in Algorithm 3. Similarly to MR4078808 , we adopt a
conservative rule of only increasing the trial stepsize $\overline{\tau}_{k}$
when two consecutive trial stepsizes were accepted without decreasing them.

1:$\gamma>1$, $\overline{\tau}_{0}>0$.

2:Obtain $\tau_{0}$ by Steps 5-7 of Algorithms 1 or 2.

3:Set $\overline{\tau}_{1}:=\max\\{\tau_{0},t_{\min}\\}$ and obtain $\tau_{1}$
by Steps 5-7 of Algorithms 1 or 2.

4:for $k=2,3,\ldots$ do

5: if $\tau_{k-2}=\overline{\tau}_{k-2}$ and
$\tau_{k-1}=\overline{\tau}_{k-1}$ then

6: $\overline{\tau}_{k}:=\gamma\tau_{k-1}$;

7: else

8: $\overline{\tau}_{k}:=\max\\{\tau_{k-1},t_{\min}\\}$.

9: end if

10: Obtain $\tau_{k}$ by Steps 5-7 of Algorithms 1 or 2.

11:end for

Algorithm 3 Self-adaptive trial stepsize

The codes in the first subsection below were written and ran in MATLAB version
R2021b, while for the second subsection we used Python 3.8. The tests were ran
on a desktop of Intel Core i7-4770 CPU 3.40GHz with 32GB RAM, under Windows 10
(64-bit).

###  6.1 Smooth DC Models in Biochemistry

Here we consider the problem motivating the development of BDCA in
AragonArtacho2018 , which consists of finding a steady state of a dynamical
equation arising in the modeling of biochemical reaction networks. We ran our
experiments on the same 14 biochemical reaction network models tested in
AragonArtacho2018 ; MR4078808 . The problem can be modeled as finding a zero
of the function

$f(x):=\left([F,R]-[R,F]\right)\exp\left(w+[F,R]^{T}x\right),$

where $F,R\in\mathbb{Z}_{\geq 0}^{m\times n}$ denote the forward and reverse
_stoichiometric matrices_ , respectively, where $w\in\mathbb{R}^{2n}$ is the
componentwise logarithm of the _kinetic parameters_ , where $\exp(\cdot)$ is
the componentwise exponential function, and where $[\,\cdot\,,\cdot\,]$ stands
for the horizontal concatenation operator. Finding a zero of $f$ is equivalent
to minimizing the function $\varphi(x):=\|f(x)\|^{2}$, which can be expressed
as a difference of the convex functions

$g(x):=2\left(\|p(x)\|^{2}+\|c(x)\|^{2}\right)\quad\text{and}\quad
h(x):=\|p(x)+c(x)\|^{2},$ (59)

where the functions $p(x)$ and $c(x)$ are given by

$p(x):=[F,R]\exp\left(w+[F,R]^{T}x\right)\quad\text{and}\quad
c(x):=[R,F]\exp\left(w+[F,R]^{T}x\right).$

In addition, it is also possible to write

$\varphi(x)=\|f(x)\|^{2}=\|p(x)-c(x)\|^{2}=\|p(x)\|^{2}+\|c(x)\|^{2}-2p(x)c(x),$

and so $\varphi(x)$ can be decomposed as the difference of the functions

$g(x):=\|p(x)\|^{2}+\|c(x)\|^{2}\quad\text{and}\quad h(x)=2p(x)c(x)$ (60)

with $g$ being convex. Therefore, $\nabla^{2}g(x)$ is $0$-lower definite, and
minimizing $\varphi$ can be tackled with Algorithm 1 by choosing
$\rho_{k}\geq\zeta$ for some fixed $\zeta>0$. As shown in AragonArtacho2018 ,
the function $\varphi$ is real analytic and thus satisfies the
Kurdyka–Łojasiewicz assumption of Theorem 4.2, but as observed in
(AragonArtacho2018, , Remark 5), a linear convergence rate cannot be
guaranteed.

Our first task in the conducted experiments was to decide how to set the
parameters $\zeta$ and $\rho_{k}$. We compared the strategy of taking
$\rho_{k}$ equal to some fixed value for all $k$, setting a decreasing
sequence bounded from below by $\zeta$, and choosing
$\rho_{k}=c\|w_{k}\|+\zeta$ for some constant $c>0$. In spite of Remark 7(ii),
$\zeta$ was added in the last strategy to guarantee both Theorem 3.1(ii) and
Theorem 4.2. We took $\zeta=10^{-8}$ and a constant $c=5$, which worked well
in all the models. We tried several options for the decreasing strategy, of
which a good choice seemed to be $\rho_{k}=\frac{\|w_{0}\|}{10^{\lfloor
k/50\rfloor}}+\zeta$, where $\lfloor\cdot\rfloor$ denotes the floor function
(i.e., the parameter was initially set to $\|w_{0}\|$ and then divided by $10$
every 50 iterations). The best option was this decreasing strategy, as can be
observed in the two models in Figure 4, and this was the choice for our
subsequent tests.

Figure 4: Comparison of the objective values for three strategies for setting
the regularization parameter $\rho_{k}$: constant (with values $10^{6}$,
$10^{5}$, $10^{3}$ and $1$), decreasing, and adaptive with respect to the
value of $\|w_{k}\|$.

###### Experiment 1

For finding a steady state of each of the 14 biochemical models, we compared
the performance of Algorithm 1 and BDCA with self-adaptive strategy, which was
the fastest method tested in MR4078808 (on average, 6.7 times faster than
DCA). For each model, 5 kinetic parameters were randomly chosen with
coordinates uniformly distributed in $(-1,1)$, and 5 random starting points
with random coordinates in $(-2,2)$ were picked. BDCA was ran using the same
parameters as in MR4078808 , while we took $\sigma=\beta=0.2$ for Algorithm 1.

We considered two strategies for setting the trial stepsize
$\overline{\tau}_{k}$ in Step 4 of Algorithm 1: constantly initially set to
50, and self-adaptive strategy (Algorithm 3) with $\gamma=2$ and
$t_{\min}=10^{-8}$. For each model and each random instance, we computed 500
iterations of BDCA with self-adaptive strategy and then ran Algorithm 1 until
the same value of the target function $\varphi$ was reached. As in
AragonArtacho2018 , the BDCA subproblems were solved by using the function
fminunc with optimoptions('fminunc', 'Algorithm', 'trust-region', 'GradObj',
'on', 'Hessian', 'on', 'Display', 'off', 'TolFun', 1e-8, 'TolX', 1e-8).

The results are summarized in Figure 5, where we plot the ratios of the
running times between BDCA with self-adaptive stepsize and Algorithm 1 with
constant trial stepsize against Algorithm 1 with self-adaptive stepsize. On
average, Algorithm 1 with self-adaptive strategy was $6.69$ times faster than
BDCA, and was $1.33$ times faster than Algorithm 1 with constant strategy. The
lowest ratio for the times of self-adaptive Algorithm 1 and BDCA was $3.17$.
Algorithm 1 with self-adaptive stepsize was only once (out of the 70
instances) slightly slower (a ratio of 0.98) than with the constant strategy.

Figure 5: Ratios of the running times of Algorithm 1 with constant stepsize
and BDCA with self-adaptive stepsize to Algorithm 1 with self-adaptive
stepsize. For each of the models, the algorithms were run using the same
random starting points. The overall average ratio is represented with a dashed
line

In Figure 6, we plot the values of the objective function for each algorithm
and also include for comparison the results for DCA and BDCA without self-
adaptive strategy. The self-adaptive strategy also accelerates the performance
of Algorithm 1. We can observe in Figure 7 that there is a correspondence
between the drops in the objective value and large increases of the stepsizes
$\tau_{k}$ (in a similar way to what was shown for BDCA in (MR4078808, , Fig.
12)).

Figure 6: Value of the objective function (with logarithmic scale) of
Algorithm 1, DCA and BDCA for two biochemical models. The value attained after
500 iterations of BDCA with self-adaptive stepsize is shown by a dashed line.

Figure 7: Comparison of the self-adaptive and the constant (with
$\overline{\tau}_{k}=50$) choices for the trial stepsizes in Step 4 of
Algorithm 1 for two biochemical models. The plots include two scales, a
logarithmic one for the objective function values and a linear one for the
stepsizes (which are represented with discontinuous lines).

###  6.2 Solving Constrained Quadratic Optimization Models

This subsection contains numerical experiments to solve problems of
constrained quadratic optimization formalized by

$\displaystyle\mbox{minimize }\;\frac{1}{2}x^{T}Qx+b^{T}x\;\text{ subject to
}\;x\in C:=\bigcup_{i=1}^{p}C_{i},$ (61)

where $Q$ is a symmetric matrix (not necessarily positive-semidefinite),
$b\in\mathbb{R}^{n}$, and $C_{1},\ldots,C_{p}\subseteq\mathbb{R}^{n}$ are
nonempty, closed, and convex sets.

When $C=\mathbb{B}_{r}(0)$ (i.e., $p=1$), this problem is referred as the
trust-region subproblem. If $Q$ is positive-semidefinite, then (61) is a
problem of convex quadratic programming. Even when $Q$ is not positive-
semidefinite, Tao and An Tao1998 showed that this particular instance of
problem (61) could be efficiently addressed with the DCA algorithm by using
the following DC decomposition:

$g(x):=\frac{1}{2}\rho\|x\|^{2}+b^{T}x+\delta_{\mathbb{B}_{r}(0)},\quad
h(x):=\frac{1}{2}x^{T}(\rho I-Q)x,$ (62)

where $\rho\geq\|Q\|_{2}$. However, this type of decomposition would not be
suitable for problem (61) when $C$ is not convex.

As shown in Subsection 5.2, problem (61) for $p\geq 1$ can be reformulated by
using FBE (51) to be tackled with Algorithm 2 with
$\lambda\in(0,\frac{1}{\|Q\|_{2}})$. Although the decomposition in (58) may
not be suitable for DCA when $Q$ is not positive-definite, it can be
regularized by adding $\frac{1}{2}\rho\|x\|^{2}$ to both $g$ and $h$ with
$\rho\geq\max\\{0,-2\lambda_{\min}(Q)\\}$. Such a regularization would
guarantee the convexity of the resulting functions $g$ and $h$ given by

$\displaystyle g(x)$
$\displaystyle:=\frac{1}{2}x^{T}\left(Q+(\rho+\lambda^{-1})I\right)x+b^{T}x,$
(63) $\displaystyle h(x)$ $\displaystyle:=\frac{1}{2}x^{T}\left(2Q+\rho
I\right)x+b^{T}x+{{\mathtt{A}}_{\lambda}\delta_{C}}((I-\lambda Q)x-\lambda
b).$ (64)

The function $g$ in (62) is not smooth, but the function $g$ in (63) is. Then
it is possible to apply BDCA MR4078808 to formulation (63)–(64) in order to
accelerate the convergence of DCA. Note that it would also be possible to do
it with (62) if the $\ell_{1}$ or $\ell_{\infty}$ balls were used; see
Artacho2019 for more details.

Let us describe two numerical experiments to solve problem (61).

###### Experiment 2

Consider (61) with $C=\mathbb{B}_{r}(0)$ and replicate the hardest setting in
Tao1998 , which was originally considered in More1983 . Specifically, in this
experiment we generated potentially difficult cases by setting $Q:=UDU^{T}$
for some diagonal matrix $D$ and orthogonal matrix $U:=U_{1}U_{2}U_{3}$ with
$U_{j}:=I-2u_{j}u_{j}^{T}/\|u_{j}\|^{2}$, $j=1,2,3$. The components of $u_{j}$
were random numbers uniformly distributed in $(-1,1)$, while the elements in
the diagonal of $D$ were random numbers in $(-5,5)$. We took $b:=Uz$ for some
vector $z$ whose elements were random numbers uniformly distributed in
$(-1,1)$ except for the component corresponding to the smallest element of
$D$, which was set to $0$. The radius $r$ was randomly chosen in the interval
$(\|d\|,2\|d\|)$, where $d_{i}:=z_{i}/(D_{ii}-\lambda_{\min}(D))$ if
$D_{ii}\neq\lambda_{\min(D)}$ and $0$ otherwise.

For each $n\in\\{100,200,\ldots,900,1000,1250,1500,\ldots,3750,4000\\}$, we
generated 10 random instances, took for each instance a random starting point
in $\mathbb{B}_{r}(0)$, and ran from it the four algorithms described above:
DCA applied to formulation (62) (without FBE), DCA and BDCA applied to
(63)–(64), and Algorithm 2. We took $\lambda=0.8/\|Q\|_{2}$ as the parameter
for FBE (both for DCA and Algorithm 2). The regularization parameter $\rho$
was chosen as $\max\\{0,-2\lambda_{\min}(Q)\\}$ for DCA with FBE and
$0.1+\max\\{0,-2\lambda_{\min}(Q)\\}$ for BDCA, as $h$ should be strongly
convex. Both Algorithm 2 and BDCA were ran with the self-adaptive trial
stepsize for the backtracking step introduced in MR4078808 with parameters
$\sigma=\beta=0.2$ and $\gamma=4$, and with $t_{\min}=10^{-6}$. For the shake
of fairness, we did not compute function values for the runs of DCA at each
iteration, since it is not required by the algorithm. Instead, we used for
both versions of DCA the stopping criterion from Tao1998 that $er\leq
10^{-4}$, where

$er=\left\\{\begin{array}[]{lc}\left\|x^{k+1}-x^{k}\right\|/\left\|x^{k}\right\|&\text{
if }\left\|x^{k}\right\|>1,\\\ \left\|x^{k+1}-x^{k}\right\|&\text{
otherwise.}\end{array}\right.$

As DCA with FBE was clearly the slowest method, we took the function value of
the solution returned by DCA without FBE as the target value for both
Algorithm 2 and BDCA, so these algorithms were stopped when that function
value was reached. In Figure 8, we plot the time ratio of each algorithm
against Algorithm 2. On average, Algorithm 2 was more than 5 times faster than
DCA with FBE and more than 2 times faster than DCA without FBE. BDCA greatly
accelerated the performance of DCA with FBE, but still Algorithm 2 was more
than 1.5 times faster. Only for size 300, the performance of DCA without FBE
was comparable to that of Algorithm 2. We observe on the right plot that the
advantage of Algorithm 2 is maintained for larger sizes.

Figure 8: Time ratio for 10 random instances of DCA with FBE, DCA without FBE,
and BDCA with respect to Algorithm 2. Average ratio within each size is
represented with a triangle for DCA with FBE, with a square for DCA without
FBE and with a circle for BDCA. The overall average ratio for each pair of
algorithms is represented by a dotted line.

###### Experiment 3

With the aim of finding the minimum of a quadratic function with integer and
box constraints, we modified the setting of Experiment 2 and considered
instead a set $C$ composed by $9^{n}$ balls of various radii centered at
$\\{-4,-3,-2,-1,0,1,2,3,4\\}^{n}$, with
$n\in\\{2,10,25,50,100,200,500,1000\\}$. As balls of radius $1/2\sqrt{n}$
cover the region $[-4,4]^{n}$, we ran our tests with balls of radii
$c/2\sqrt{n}$ with $c\in\\{0.1,0.2,\ldots,0.8,0.9\\}$. This time we considered
both convex and nonconvex objective functions. The nonconvex case was
generated as in Experiment 2, while for the convex case, the elements of the
diagonal of $D$ were chosen as random numbers uniformly distributed in
$(0,5)$.

For each $n$ and $r$, 100 random instances were generated. For each instance,
a starting point was chosen with random coordinates uniformly distributed in
$[-5,5]^{n}$. As the constraint sets are nonconvex, FBE was also needed to run
DCA. The results are summarized in Table 3, where for each $n$ and each
radius, we counted the number of instances (out of 100) in which the value of
$\varphi_{\lambda}$ at the rounded output of DCA and BDCA was lower and higher
than that of Algorithm 2 when ran from the same starting point. We used the
same parameter settings for the algorithms as in Experiment 2. Finally, we
plot in Figure 9 two instances in $\mathbb{R}^{2}$ in which Algorithm 2
reached a better solution.

| Radius of the balls
---|---
7pt. |  Alg. 2 vs | $\frac{1}{20}\sqrt{n}$ | $\frac{2}{20}\sqrt{n}$ | $\frac{3}{20}\sqrt{n}$ | $\frac{4}{20}\sqrt{n}$ | $\frac{5}{20}\sqrt{n}$ | $\frac{6}{20}\sqrt{n}$ | $\frac{7}{20}\sqrt{n}$ | $\frac{8}{20}\sqrt{n}$ | $\frac{9}{20}\sqrt{n}$
$n=2$ | DCA | 0/34 | 1/20 | 1/26 | 1/19 | 0/19 | 0/18 | 0/3 | 1/1 | 0/2
BDCA | 6/12 | 2/13 | 3/14 | 4/9 | 0/4 | 1/10 | 0/2 | 1/1 | 0/2
$n=10$ | DCA | 2/89 | 2/83 | 1/66 | 5/53 | 8/28 | 3/7 | 1/1 | 3/1 | 1/0
BDCA | 21/68 | 38/53 | 33/39 | 23/24 | 18/16 | 2/7 | 0/0 | 0/1 | 0/0
$n=25$ | DCA | 0/99 | 0/98 | 2/87 | 11/58 | 9/32 | 3/8 | 2/9 | 2/2 | 5/4
BDCA | 16/83 | 29/71 | 40/58 | 37/40 | 13/26 | 2/3 | 0/1 | 0/0 | 1/2
$n=50$ | DCA | 0/100 | 0/100 | 0/91 | 2/86 | 13/41 | 14/12 | 9/12 | 6/10 | 12/12
BDCA | 8/92 | 6/94 | 31/69 | 36/53 | 16/28 | 8/8 | 6/5 | 3/4 | 5/3
$n=100$ | DCA | 0/100 | 0/100 | 0/99 | 9/87 | 18/49 | 18/31 | 12/22 | 18/20 | 11/21
BDCA | 2/98 | 6/94 | 39/61 | 36/61 | 23/33 | 16/14 | 9/8 | 9/8 | 13/9
$n=200$ | DCA | 0/100 | 0/100 | 0/100 | 1/98 | 23/64 | 31/41 | 25/29 | 22/30 | 20/41
BDCA | 3/97 | 2/98 | 38/62 | 37/63 | 33/39 | 27/17 | 18/18 | 14/13 | 16/18
$n=500$ | DCA | 0/100 | 0/100 | 0/100 | 1/99 | 6/94 | 15/80 | 27/61 | 29/65 | 36/48
BDCA | 0/100 | 1/99 | 41/59 | 44/56 | 33/63 | 25/56 | 34/39 | 32/47 | 17/35

(a) Convex case

| Radius of the balls
---|---
7pt. |  Alg. 2 vs | $\frac{1}{20}\sqrt{n}$ | $\frac{2}{20}\sqrt{n}$ | $\frac{3}{20}\sqrt{n}$ | $\frac{4}{20}\sqrt{n}$ | $\frac{5}{20}\sqrt{n}$ | $\frac{6}{20}\sqrt{n}$ | $\frac{7}{20}\sqrt{n}$ | $\frac{8}{20}\sqrt{n}$ | $\frac{9}{20}\sqrt{n}$
$n=2$ | DCA | 1/8 | 1/9 | 1/10 | 1/6 | 0/6 | 0/6 | 0/8 | 3/2 | 1/0
BDCA | 2/4 | 1/4 | 1/3 | 3/4 | 0/5 | 0/4 | 0/6 | 3/2 | 1/0
$n=10$ | DCA | 9/39 | 4/39 | 7/39 | 4/35 | 10/30 | 3/27 | 5/45 | 2/34 | 8/29
BDCA | 9/31 | 11/33 | 13/29 | 6/31 | 11/29 | 6/25 | 7/38 | 5/29 | 10/30
$n=25$ | DCA | 6/69 | 13/67 | 7/62 | 5/61 | 10/53 | 3/59 | 6/56 | 3/72 | 3/66
BDCA | 16/58 | 16/63 | 16/55 | 12/48 | 9/52 | 11/52 | 13/52 | 12/57 | 11/58
$n=50$ | DCA | 11/81 | 10/79 | 8/87 | 5/90 | 3/87 | 4/80 | 2/86 | 5/89 | 8/81
BDCA | 24/68 | 21/64 | 23/70 | 17/73 | 14/75 | 9/73 | 10/75 | 18/74 | 16/71
$n=100$ | DCA | 4/96 | 6/94 | 4/94 | 5/94 | 4/96 | 3/97 | 2/98 | 7/91 | 9/91
BDCA | 15/85 | 16/83 | 18/80 | 14/84 | 17/83 | 11/89 | 9/91 | 20/79 | 19/80
$n=200$ | DCA | 4/96 | 4/96 | 4/96 | 2/98 | 1/99 | 2/98 | 4/96 | 3/97 | 0/100
BDCA | 11/89 | 16/84 | 11/89 | 8/92 | 6/94 | 11/89 | 10/90 | 13/87 | 8/92
$n=500$ | DCA | 1/99 | 2/98 | 0/100 | 0/100 | 0/100 | 1/99 | 1/99 | 2/98 | 1/99
BDCA | 12/88 | 17/83 | 15/85 | 9/91 | 15/85 | 11/89 | 9/91 | 18/82 | 20/80

(b) Nonconvex case (c) For different values of $n$ (space dimension) we
computed 100 random instances of problem (61) with $Q$ positive definite and
$C$ formed by the union of balls whose centers have integer coordinates
between $-4$ and $4$. We counted the number of instances in which DCA and BDCA
obtained a lower/upper value than Algorithm 2.

Figure 9: Two instances of problem (61). On the left, both line searches of
Algorithm 2 and BDCA help to reach a better solution for a nonconvex quadratic
function, while only Algorithm 2 succeeds on the right for the convex case.

##  7 Conclusion and Future Research

This paper proposes and develops a novel RCSN method to solve problems of
difference programming whose objectives are represented as differences of
generally nonconvex functions. We establish well-posedness of the proposed
algorithm and its global convergence under appropriate assumptions. The
obtained results exhibit advantages of our algorithm over known algorithms for
DC programming when both functions in the difference representations are
convex. We also develop specifications of the main algorithm in the case of
structured problems of constrained optimization and conduct numerical
experiments to confirm the efficiency of our algorithms in solving practical
models.

In the future research, we plan to relax assumptions on the program data
ensuring the linear, superlinear, and quadratic convergence rates for RCSN and
also extend the spectrum of applications to particular classes of constrained
optimization problems as well as to practical modeling.

## References

  * (1) Aragón-Artacho, F.J., Goberna, M.A., López, M.A., Rodríguez, M.M.L.: Nonlinear optimization. Springer, Cham (2019)
  * (2) Aragón-Artacho, F.J., Geoffroy, M.H.: Metric subregularity of the convex subdifferential in Banach spaces. J. Nonlinear Convex Anal. 15, 35–47 (1014)
  * (3) Aragón-Artacho, F.J., Campoy, R., Vuong, P.T.: Using positive spanning sets to achieve d-stationarity with the boosted DC algorithm. Vietnam J. Math. 48, 363–376 (2020)
  * (4) Aragón-Artacho, F.J., Campoy, R., Vuong, P.T.: The boosted DC algorithm for linearly constrained DC programming. Set-Valued Var. Anal. 30, 1265–1289 (2022)
  * (5) Aragón-Artacho, F.J., Fleming, R.M.T., Vuong, P.T.: Accelerating the DC algorithm for smooth functions. Math. Program. 169, 95–118 (2018)
  * (6) Aragón-Artacho, F.J., Vuong, P.T.: The boosted difference of convex functions algorithm for nonsmooth functions. SIAM J. Optim. 30, 980–1006 (2020)
  * (7) Asplund, E.: Fréchet differentiability of convex functions. Acta Math. 121, 31–47 (1968).
  * (8) Bernard, F., Thibault, L.: Prox-regularity of functions and sets in Banach spaces. Set-Valued Anal. 12, 25–47 (2004)
  * (9) Bernard, F., Thibault, L.: Uniform prox-regularity of functions and epigraphs in Hilbert spaces. Nonlinear Anal. 60, 187–207 (2005)
  * (10) Colombo, G., Henrion, R., Hoang, N.D.; Mordukhovich, B.S.: Optimal control of sweeping processes over polyhedral control sets. J. Diff. Eqs. 260, 3397–3447 (2016)
  * (11) de Oliveira, W.: The ABC of DC programming. Set-Valued Var. Anal. 28, 679–706 (2020)
  * (12) Ding, C., Sun, D., Ye, J.J.: First-order optimality conditions for mathematical programs with semidefinite cone complementarity constraints. Math. Program. 147, 539–379 (2014)
  * (13) Drusvyatskiy, D., Mordukhovich, B.S., Nghia, T.T.A.: Second-order growth, tilt stability, and metric regularity of the subdifferential. J. Convex Anal. 21, 1165–1192 (2014)
  * (14) Facchinei, F., Pang, J.-S.: Finite-Dimensional Variational Inequalities and Complementarity Problems, I, II. Springer, New York (2003)
  * (15) Gfrerer, H., Outrata, J.V.: On a semismooth∗ Newton method for solving generalized equations. SIAM J. Optim. 31, 489–517 (2021)
  * (16) Henrion, R., Mordukhovich, B.S., Nam, N.M.: Second-order analysis of polyhedral systems in finite and infinite dimensions with applications to robust stability of variational inequalities. SIAM J. Optim. 20, 2199–2227 (2010)
  * (17) Henrion, R., Outrata, J., Surowiec, T.: On the co-derivative of normal cone mappings to inequality systems. Nonlinear Anal. 71, 1213–1226 (2009)
  * (18) Henrion, R., Römisch, W.: On $M$-stationary points for a stochastic equilibrium problem under equilibrium constraints in electricity spot market modeling. Appl. Math. 52, 473–494 (2007)
  * (19) Hiriart-Urruty, J.-B.: Generalized differentiability, duality and optimization for problems dealing with differences of convex functions. In: Ponstein, J. (ed.) Convexity and Duality in Optimization. Lecture Notes Econ. Math. Syst. 256, pp. 37–70. Springer, Berlin (1985)
  * (20) Izmailov, A.F., Solodov, M.V.: Newton-Type Methods for Optimization and Variational Problems. Springer, Cham (2014)
  * (21) Khanh, P.D., Mordukhovich, B.S., Phat, V.T.: A generalized Newton method for subgradient systems. Math. Oper. Res. (2022), DOI 10.1287/moor.2022.1320
  * (22) Khanh, P.D., Mordukhovich, B.S., Phat, V.T., Tran, D.B.: Generalized Newton algorithms in nonsmooth optimization via second-order subdifferentials. J. Global Optim. (2022). DOI 10.1007/s10898-022-01248-7
  * (23) Khanh, P.D., Mordukhovich, B.S., Phat, V.T., Tran, D.B.: Globally convergent coderivative-based generalized Newton methods in nonsmooth optimization (2022). arXiv:2109.02093
  * (24) Li, W., Bian, W., Toh, K.-C.: Difference-of-convex algorithms for a class of sparse group $\ell_{0}$ regularized optimization problems. SIAM J. Optim. 32, 1614–1641 (2022)
  * (25) Mordukhovich, B.S.: Sensitivity analysis in nonsmooth optimization. In: Field, D.A., Komkov, V.(eds) Theoretical Aspects of Industrial Design, pp. 32–46. SIAM Proc. Appl. Math. 58. Philadelphia, PA (1992)
  * (26) Mordukhovich, B.S.: Variational Analysis and Generalized Differentiation, I: Basic Theory, II: Applications. Springer, Berlin (2006)
  * (27) Mordukhovich, B.S.: Variational Analysis and Applications. Springer, Cham (2018)
  * (28) Mordukhovich, B.S., Outrata, J.V.: On second-order subdifferentials and their applications. SIAM J. Optim. 12, 139–169 (2001)
  * (29) Mordukhovich, B.S., Rockafellar, R.T.: Second-order subdifferential calculus with applications to tilt stability in optimization. SIAM J. Optim. 22, 953–986 (2012)
  * (30) Mordukhovich, B.S., Sarabi, M.E.: Generalized Newton algorithms for tilt-stable minimizers in nonsmooth optimization. SIAM J. Optim. 31, 1184–1214 (2021)
  * (31) Moré, J.J., Sorensen, D.C.: Computing a trust region step. SIAM J. Sci. Statist. Comput. 4. 553–572 (1983)
  * (32) Ostrowski, A.M.: Solution of Equations and Systems of Equations, 2nd ed. Academic Press, Cambridge, MA (1966)
  * (33) Outrata, J.V., Sun, D.: On the coderivative of the projection operator onto the second-order cone. Set-Valued Anal. 16 (999-1014 (2008)
  * (34) Patrinos, P., Bemporad, A.: Proximal Newton methods for convex composite optimization. In: 52nd IEEE Conf. Dec. Cont., pp. 2358–2363. Florence, Italy (2013)
  * (35) Rockafellar, R.T., Wets, R.J-B.: Variational Analysis. Springer, Berlin (1998)
  * (36) Tao, P.D., An, L.T.H.: Convex analysis approach to DC programming: theory, algorithms and applications. Acta Math. Vietnam. 22, 289–355 (1997)
  * (37) Tao, P.D., An, L.T.H.: A DC optimization algorithm for solving the trust-region subproblem. SIAM J. Optim. 8, 476–505 (1998)
  * (38) Tao, P.D., Bernoussi, E.S.: Algorithms for solving a class of nonconvex optimization problems. Methods of subgradients. North-Holland Math. Stud. 129, 249–271 (1986)
  * (39) Themelis A., Stella, L. and Patrinos, P.: Forward-backward envelope for the sum of two nonconvex functions: further properties and nonmonotone linesearch algorithms. SIAM J. Optim. 28, 2274–2303 (2018)
  * (40) Toland, J.F.: On subdifferential calculus and duality in non-convex optimization. Mem. Soc. Math. France. 60, 177–183 (1979)
  * (41) Truong, T.T., Nguyen, H.T.: Backtracking gradient descent method and some applications in large scale optimisation, II: Algorithms and experiments. Appl. Math. Optim. 84, 2557–2586 (2021)
  * (42) Yao, J.-C., Yen, N.D.: Coderivative calculation related to a parametric affine variational inequality. Part 1: Basic calculation. Acta Math. Vietnam. 34, 157–172 (2009)