∎

11institutetext: $1$ Mines Paris, Université PSL, CMA – Centre de
Mathématiques Appliquées, France
11email<EMAIL_ADDRESS>
$2$ University of Southern California

# An implementable proximal-type method for computing critical points to
minimization problems with a nonsmooth and nonconvex constraint

Gregorio M. Sempere1 Welington de Oliveira1 Johannes O. Royset2

###### Abstract

This work proposes an implementable proximal-type method for a broad class of
optimization problems involving nonsmooth and nonconvex objective and
constraint functions. In contrast to existing methods that rely on an ad hoc
model approximating the nonconvex functions, our approach can work with a
nonconvex model constructed by the pointwise minimum of finitely many convex
models. The latter can be chosen with reasonable flexibility to better fit the
underlying functions’ structure. We provide a unifying framework and analysis
covering several subclasses of composite optimization problems and show that
our method computes points satisfying certain necessary optimality conditions,
which we will call model criticality. Depending on the specific model being
used, our general concept of criticality boils down to standard necessary
optimality conditions. Numerical experiments on some stochastic reliability-
based optimization problems illustrate the practical performance of the
method.

###### Keywords:

Composite Optimization Nonsmooth Optimization Nonconvex Optimization
Variational Analysis

###### MSC:

49J52 49J53 49K99 90C26

††journal: -

##  1 Introduction

Modern optimization problems, such as those in stochastic optimization,
machine learning, uncertainty quantification, and others, often involve
functions that fail to be smooth and convex. Nevertheless, many such problems
possess some structure that can be exploited in numerical optimization.
Existing nonsmooth and nonconvex optimization methods typically rely on an ad
hoc model exploiting a specific (composite) structure Burke_1995 ;
Sagastizabal_2013 ; Lewis_2016 ; Cui_Pang_Sen_2018 ; Liu_Cui_2020 ;
Cui_Pang_Liu_2022 ; AragonPerezTorregrosa_2023 . Most such methods deal with
the convex-constrained setting. For instance, Burke_1985 ; Burke_1995 ;
Sagastizabal_2012  exploits convexity and smoothness of certain components
functions yielding the composite one, Lewis_2016  and
AragonPerezTorregrosa_2023  profit from a certain degree of smoothness,
Pham_2005 ; Oliveira_2019 ; Cui_Pang_Sen_2018  from a Difference-of-Convex
structure, and Ackooij_Oliveira_2020  exploits the difference-of-locally-
Lipschitz format. A framework for constructing and analyzing approximations of
optimization problems with a composite structure is given in Royset_2023 , and
(Cui_Pang_2022, , Chapter 7) presents several algorithms based on the
principle of surrogation.

This work considers a broader class of composite optimization problems
involving nonsmooth and nonconvex objective and constraint functions. The
investigated setting exploits the available structure so that the problem can
be locally approximated by a model even if subgradients are unavailable. This
is the case for some composite functions whose first-order information is
available for the components but not necessarily for the composition.
Applications include Difference-of-Convex (DC)-constrained DC problems,
classification, and sample average approximations of chance-constrained and
buffered probability optimization problems.

Some optimization algorithms exist for handling different families of problems
fitting into our setting. For instance, Apkarian2008  deals with the so-called
mixed $H_{2}/H_{\infty}$ control problem, Cui_Pang_Liu_2022  focuses on the
particular case of chance-constrained problems, and the DC-constrained case is
investigated in Pham_Penalization_2012 ; Strekalovsky_Minarchenko_2017 ;
Pang_2017 ; Montonen_Joki_2018 ; Cui_2018 ; Javal_2021 . A broader setting
than DC was recently investigated in Ksenia_2023 , where the objective and
constraint functions can be expressed as the difference of convex and weakly-
convex functions. These papers build different models to provide implementable
algorithms in their respective settings. The main idea is to exploit the
functions’ structure to approximate the problem locally. For instance,
complicating (component) functions are linearized in Pham_Penalization_2012 ;
Strekalovsky_Minarchenko_2017 ; Pang_2017 ; Montonen_Joki_2018 ; Cui_2018 ;
Javal_2021 ; Ksenia_2023  whereas higher-order approximations are considered
in AragonPerezMordukhovich_2023  under structured nonsmoothness. In addition,
these models must accommodate nonconvex constraints. Standard ways to face
constrained optimization problems are through penalization
Pham_Penalization_2012 ; Byun_Oliveira_Royset_2023  and (Cui_Pang_2022, , §9),
linearizations of constraints Pang_2017 ; Cui_2018 ; Ackooij_Oliveira_OMS_2019
, and, as considered in this work, the so-called improvement function
Montonen_Joki_2018 ; Javal_2021 ; Sagastizabal_Solodov_2005 ; Apkarian2008
that jointly combines objective and constraint functions.

Each of these references develops proximal-type algorithms to compute a point
satisfying some necessary optimality condition denoted by criticality. Many
definitions of criticality exist in nonsmooth and nonconvex optimization,
which are closely related to the working model for the problem. Therefore,
each of these references exploits a model and provides a dedicated algorithm
and convergence analysis. Differently, this work proposes a general yet
implementable proximal-type algorithm with guaranteed convergence to critical
points. Our study is comprehensive, and the type of criticality we get depends
explicitly on the model’s properties. The latter can be chosen with reasonable
freedom as long as some conditions hold to ensure its tractability.

In a nutshell, our approach is based on a nonconvex model (for the improvement
function) constructed by the pointwise minimum of finitely many _convex
models_ , which are obtained using the available information about the
functions. For instance, if some of the available oracles for the component
functions can provide us with more than one subgradient for every given point,
we might benefit from this fact to construct more than one convex model
approximating the functions locally. The pointwise minimum of such convex
models leads to nonconvex but computationally solvable subproblems yielding
new iterates. Thanks to this more general model, our method and analysis cover
several classes of (nonlinearly constrained) composite optimization problems.
Our analysis is based on the subdifferentiability (see (Clarke_1987, , §2),
(RoysetWets_2022, , §4) and (Cui_Pang_2022, , §4)) and epigraphical nesting
(see (RoysetWets_2022, , §6.C) and HigleSen_1992 ) theories. The latter is a
known relaxation of epi-convergence (see (RoysetWets_2022, , §5) and
(Rockafellar_Wets_2009, , §7)) with a broader field of use that allows us to
accommodate the convergence analysis of implementable algorithms, such the one
proposed in Section 4 below.

The manuscript is organized as follows. Section 2 recalls standard definitions
and fundamental concepts for our development, such as the improvement function
and epigraphical nesting. Section 3 presents the basic assumptions of the
working model, together with the type of criticality to attain. Our algorithm
is given in Section 4 and its analysis relies on epigraphical nesting.
Sections 5, 6, and 7 present some examples of classes of functions for which
we can provide implementable models for our algorithm. The manuscript closes
with Section 8 illustrating the numerical performance of our approach with two
stochastic optimization problems.

#### Notation:

We consider $f:\mathbb{R}^{n}\to\mathbb{R}$ as a real-valued function, while
we use $\mathsf{f}:\mathbb{R}^{n}\to\mathbb{R}\cup\\{+\infty\\}$ to refer to
an extended real-valued function with a nonempty domain. We denote by $X$ a
convex and closed set and by ${\rm i}_{X}$ the indicator function of $X$, that
equals $0$ if $x\in X$ and equals $+\infty$ otherwise. Also, we denote by
${\boldsymbol{1}}_{A}$ the characteristic function of $A$, that equals $1$ if
$x\in A$ and equals $0$ otherwise. The normal cone of the _(closed and
convex)_ set $X$ at $\bar{x}\in X$ is

$N_{X}(\bar{x})=\\{\xi\in\mathbb{R}^{n}:\langle\xi,x-\bar{x}\rangle\leq
0,\;\forall x\in X\\},$

and $N_{X}(\bar{x})=\emptyset$ if $\bar{x}\notin X$. Given an optimization
problem such as $\min_{x\in X}f(x)$, where $f$ is locally Lipschitz, we say
that $\bar{x}\in X$ is a stationary point if $0\in\partial^{\tt
C}f(\bar{x})+N_{X}(\bar{x})$, where

$\partial^{\tt C}f(\bar{x})={\tt
co}\\{\xi\in\mathbb{R}^{n}:\exists\;x^{k}\to\bar{x},\nabla f(x^{k})\to\xi\\}$

is the Clarke subdifferential of $f$ at $\bar{x}$. If $f$ happens to be
convex, it is well-known that $\partial^{\tt C}f(\bar{x})=\partial
f(\bar{x})$, where

$\partial f(\bar{x})=\\{\xi\in\mathbb{R}^{n}:\langle\xi,x-\bar{x}\rangle\leq
f(x)-f(\bar{x}),\forall x\in\mathbb{R}^{n}\\}.$

We consider the tangent cone to a closed set $A$ at a point $\bar{x}\in A$ in
the Bouligand sense:

$\mathcal{T}_{A}(\bar{x})=\left\\{d\in\mathbb{R}^{n}:\exists\\{y^{k}\\}_{k}\subseteq
A\;,t_{k}>0\;,y^{k}\to\bar{x}\;,\;t_{k}\to
0:d=\lim_{k\to\infty}\frac{y^{k}-\bar{x}}{t_{k}}\right\\}.$

We denote the epigraph of a function
$\mathsf{f}:\mathbb{R}^{n}\to\mathbb{R}\cup\\{+\infty\\}$ by
$\operatorname{epi}\left(\mathsf{f}\right)=\\{(x,\alpha)\in\mathbb{R}^{n+1}:\mathsf{f}(x)\leq\alpha\\}$,
and the inner and outer limits of sequence of sets
$D_{k}\subseteq\mathbb{R}^{n}$ by
$\operatorname{LimInn}\left(D_{k}\right)=\\{\bar{x}\in\mathbb{R}^{n}:\exists
y^{k}\in D_{k},y^{k}\to\bar{x}\\}$ and
$\operatorname{LimOut}\left(D_{k}\right)=\\{\bar{x}\in\mathbb{R}^{n}:\exists\mathcal{K}\subset\mathbb{N},\;y^{k}\in
D_{k},\lim_{k\in\mathcal{K}\to\infty}y^{k}=\bar{x}\\}$.

##  2 Problem statement and prerequisites

Given locally Lipschitz continuous functions $f,c:\mathbb{R}^{n}\to\mathbb{R}$
and $\emptyset\neq X\subseteq\mathbb{R}^{n}$, a convex and closed set, we are
interested in optimization problems such as

$\min_{x\in X}\;f(x)\quad\text{s.t.}\quad c(x)\leq 0.$ (1)

We assume that the domains of both $f$ and $c$ are the whole space
$\mathbb{R}^{n}$, however, this could be relaxed to the set $X$ being
contained in the interior of both domains. Functions $f$ and $c$ need not be
convex or differentiable, which makes this problem challenging from the
numerical optimization point of view unless some favorable structure is
assumed. Even when the problem is well-structured, e.g., DC problems,
computing a solution with a certificate of global optimality is out of reach
in real-life applications such as those coming from stochastic optimization or
machine learning. This is why much research focuses on local-solution methods,
i.e., algorithms capable of computing a point satisfying certain necessary
optimality conditions, denoted by criticality. What we mean by criticality is
specified in Section 3 below.

####  Improvement function.

We now recall an important mathematical tool in nonlinearly-constrained
nonsmooth optimization, the so-called improvement function
Sagastizabal_Solodov_2005 :

$H(y;x)=\max\\{f(y)-\tau_{f}(x)\;,\;c(y)\\}\quad\mbox{with}\quad\tau_{f}(x)=f(x)+\rho[c(x)]_{+}$
(2)

and given parameter $\rho\geq 0$. We highlight that dealing with a difficult
constraint via an improvement function is a well-known and successful strategy
in the nonsmooth optimization literature Sagastizabal_Solodov_2005 ;
Apkarian2008 ; Montonen_Joki_2018 ; Javal_2021 ; KseniaPaul_2023 . It is an
important alternative to penalty approaches, which generally require laborious
work of (penalty) parameter tuning. In contrast, the parameter $\rho$ is
simple to tune as it can be chosen such that the magnitudes of both functions
$f$ and $c$ are close111 For instance, a possible choice for $\rho$ is
$\rho=|f(x^{0})|/(1+|c(x^{0})|)$, with $x^{0}\in X$ being a given initial
point to the optimization algorithm.. The following properties of $H$, which
are trivial extensions of (Apkarian2008, , Lem. 5.1) to our setting, are
handy.

######  Lemma 1 (Lem. 5.1 in Apkarian2008  and Lem. 1 in Javal_2021 .)

The following hold.

  1. (i)

Let $\mu\geq 0$ and $x,y$ be two points in $X$ satisfying
$H(y;x)+\frac{\mu}{2}\textstyle{\|}y-x\textstyle{\|}^{2}\leq H(x;x)$:

    1. (i.1)

If $c(x)\leq 0$, then
$c(y)+\frac{\mu}{2}\textstyle{\|}y-x\textstyle{\|}^{2}\leq 0$ and
$f(y)+\frac{\mu}{2}\textstyle{\|}y-x\textstyle{\|}^{2}\leq f(x)$.

    2. (i.2)

If $c(x)>0$, then $c(y)+\frac{\mu}{2}\textstyle{\|}y-x\textstyle{\|}^{2}\leq
c(x)$.

  2. (ii)

If $\bar{x}\in X$ is a local solution to (1), then $\bar{x}$ is also a local
minimizer of $H(x;\bar{x})$ over $X$.

$\blacktriangleleft$

####  Epi-convergence and epigraphical nesting.

We now recall the concept of _epigraphical nesting_ , which is a relaxation of
_epi-convergence_ (Rockafellar_Wets_2009, , §7). Epigraphical nesting was
proposed in HigleSen_1992  as a tool for analyzing minimization algorithms
(see also (RoysetWets_2022, , §6.C)).

######  Definition 1 (Def. 7.1 in Rockafellar_Wets_2009  and Thm. 4.15 in
RoysetWets_2022 )

Given
$\mathsf{f}^{k},\mathsf{f}:\mathbb{R}^{n}\to\mathbb{R}\cup\\{+\infty\\}$, the
sequence $\\{\mathsf{f}^{k}\\}_{k}$ is said to epi-converge to $\mathsf{f}$,
denoted by $\mathsf{f}^{k}\xrightarrow{e}\mathsf{f}$, if the following
condition holds:

$\operatorname{LimOut}\left(\operatorname{epi}\left(\mathsf{f}^{k}\right)\right)\subseteq\operatorname{epi}\left(\mathsf{f}\right)\subseteq\operatorname{LimInn}\left(\operatorname{epi}\left(\mathsf{f}^{k}\right)\right).$

In that case, we refer to $\mathsf{f}$ as the epi-limit of the sequence
$\\{\mathsf{f}^{k}\\}_{k}$.$\blacktriangleleft$

Note that, by definition of inner and outer limits, the above inclusions are
equivalent to all three sets to be equal. We choose the above characterization
of epi-convergence as it gives a clear view of how it relates to epigraphical
nesting. However, expressing these inclusions in terms of convergence of
$\\{\mathsf{f}^{k}\\}_{k}$ is handy for the following proofs.

###### Lemma 2

Given $f^{k},f:\mathbb{R}^{n}\to\mathbb{R}$ and $X\subseteq\mathbb{R}^{n}$ a
closed set, it holds:

  1. (i)

$\operatorname{LimOut}\left(\operatorname{epi}\left(f^{k}+{\rm
i}_{X}\right)\right)\subseteq\operatorname{epi}\left(f+{\rm i}_{X}\right)$ if,
and only if, for all $x\in X$ and for each sequence $\\{x^{k}\\}_{k}\subseteq
X$ such that $x^{k}\to x$,

$f(x)\leq\liminf_{k}f^{k}(x^{k}).$

  2. (ii)

$\operatorname{epi}\left(f+{\rm
i}_{X}\right)\subseteq\operatorname{LimInn}\left(\operatorname{epi}\left(f^{k}+{\rm
i}_{X}\right)\right)$ if, and only if, for all $x\in X$ there exists a
sequence $\\{y^{k}\\}_{k}\subseteq X$ such that $y^{k}\to x$ and

$\limsup_{k}f^{k}(y^{k})\leq f(x).$

Proof Let us recall (Rockafellar_Wets_2009, , Prop. 7.2), which states
essentially the same result for extended-valued functions. In that setting,
there is no need to focus on points and subsequences in $X$. When applying
that result to our case, it turns out that we do not need to investigate
$x\notin X$, as both inequalities

$\limsup_{k}f^{k}(x^{k})+{\rm i}_{X}(x^{k})\leq f(x)+{\rm
i}_{X}(x)\leq\liminf_{k}f^{k}(x^{k})+{\rm i}_{X}(x^{k})$

would hold for any sequence $x^{k}\to x$, as ${\rm i}_{X}(x^{k})={\rm
i}_{X}(x)=+\infty$ eventually. Let us now assume that $x\in X$. To prove item
(i), note that $\operatorname{LimOut}\left(\operatorname{epi}\left(f^{k}+{\rm
i}_{X}\right)\right)\subseteq\operatorname{epi}\left(f+{\rm i}_{X}\right)$ is
equivalent to $\liminf_{k}f^{k}(x^{k})+{\rm i}_{X}(x^{k})\geq f(x)+{\rm
i}_{X}(x)$ for all $x^{k}\to x$ (due to (Rockafellar_Wets_2009, , Prop. 7.2)),
if there were a subsequence that happened to be out of $X$, the inequality
would hold trivially. Thus, it is enough to consider only for those sequences
$\\{x^{k}\\}_{k}\subseteq X$.

Consider now item (ii). It is stated in (Rockafellar_Wets_2009, , Prop. 7.2)
that $\operatorname{epi}\left(f+{\rm
i}_{X}\right)\subseteq\operatorname{LimInn}\left(\operatorname{epi}\left(f^{k}+{\rm
i}_{X}\right)\right)$ is equivalent to $\limsup_{k}f^{k}(y^{k})+{\rm
i}_{X}(y^{k})\leq f(x)+{\rm i}_{X}(x)$ for a given $y^{k}\to x$. If the tail
of $\\{y^{k}\\}_{k}$ were not eventually contained in $X$, the superior limit
would be infinity, violating the required inequality. Thus, the sequence
$\\{y^{k}\\}_{k}\subseteq X$ without loss of generality, which in turn proves
the statement in item (ii). $\square$

Note that for the sequence $\\{y^{k}\\}_{k}$ above, we have
$\lim_{k}f^{k}(y^{k})=f(x)$. Useful properties can be proven for a sequence of
functions that is epi-convergent (see (Rockafellar_Wets_2009, , §7)). For
instance, all cluster points of the sequence of approximate minimizers of
$f^{k}+{\rm i}_{X}$ minimize $f+{\rm i}_{X}$ under mild assumptions
(Rockafellar_Wets_1986, , Thm. 7.31 (b)). Although the importance of this
result in optimization, algorithmically designing an epi-convergent sequence
of functions is not always a simple task. For instance, it is well known that
even basic convergent algorithms, such as the cutting-plane method in convex
optimization, fail to construct a sequence of functions that epi-converges to
the objective function. The main obstacle in the concept of epi-convergence is
that it requires a sequence of functions for which we have both over- and
under-estimator subsequences around every point. It was shown in HigleSen_1992
that, for algorithmic design purposes, it is enough to have only the right-
most inclusion in Definition 1, i.e.,

$\operatorname{epi}\left(\mathsf{f}\right)\subseteq\operatorname{LimInn}\left(\operatorname{epi}\left(\mathsf{f}^{k}\right)\right).$

This condition is known as _epigraphical nesting_ , and it ensures that the
sequence of functions estimates $\mathsf{f}$ from below. This concept is
clearly weaker than epi-convergence, but it still proves useful in
optimization (see (RoysetWets_2022, , Prop. 6.14)).

##  3 Model and necessary optimality conditions

In nonconvex nonsmooth optimization, many definitions of stationary points
exist. A well-known definition is Bouligand stationarity Pang_2007 : a point
$\bar{x}\in X^{c}:=\\{x\in X:\;c(x)\leq 0\\}$ is $B$-stationary to problem (1)
if the directional derivative of $f$ with respect to $d$ is non-negative for
all $d\in\mathcal{T}_{X^{c}}(\bar{x})$. Obtaining an explicit expression of
$\mathcal{T}_{X^{c}}(\bar{x})$ is not simple, let alone an implementable one.
We refer the interested reader to Pang_2017  and Ackooij_Oliveira_2022  for
implementable expressions of $\mathcal{T}_{X^{c}}(\bar{x})$ in particular
settings of DC programming, see also Section 6 below. For this reason,
different necessary optimality conditions must come into play.

Given item (ii) of Lemma 1, we could then use
$\bar{x}\in\operatorname{argmin}_{x\in X}H(x;\bar{x})$ as a necessary
optimality condition for $\bar{x}\in X$ to solve problem (1). However, for the
class of problems we have in mind, working directly with the improvement
function is still challenging. The reason is that it may happen in practice
that (first-order) oracles for $f$ and $c$ are not available, but for
functions composing them; see Sections 5-7 below. Our main assumption is that
we can construct, with the available oracle information, a model
$\mathcal{M}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ for $H$
satisfying the following hypothesis.

######  Assumption 1 (Model’s requirements.)

For all $x\in X$ fixed, let $A(x)$ be a finite index set and
$\mathcal{M}_{a}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ be a real-
valued function such that $\mathcal{M}_{a}(\cdot;x)$ is convex for all $a\in
A(x)$. We take $\mathcal{M}(\cdot;x):=\min_{a\in
A(x)}\;\mathcal{M}_{a}(\cdot;x)$ as a model for $H(\cdot;x)$ in (2) provided
that the following conditions hold:

  1. (a)

$\mathcal{M}(x;x)=H(x;x)$ for all $x\in X$;

  2. (b)

there exists a constant $\bar{\mu}\in[0,+\infty)$ such that
$H(y;x)\leq\mathcal{M}(y;x)+\frac{\bar{\mu}}{2}\textstyle{\|}y-x\textstyle{\|}^{2}$
for all $y,x\in X$.

$\blacktriangleleft$

The hypotheses on $A(x)$ and $\mathcal{M}_{a}$ are natural and helpful when
$H$ can be written as the pointwise minimum of finitely many functions or when
oracles assessing component functions in $H$ can provide more than one
subgradient. If none are the case, we can take $A(x)$ being a singleton so
$\mathcal{M}(\cdot;x)$ is a convex model. The case $|A(x)|>1$ gives rise to a
nonconvex (DC) model. Assumption (a) is also a natural one and (b) requests
the model, plus a quadratic term, to overestimate the improvement function.
Although this last condition looks too restrictive, the following sections
show otherwise: the modulus $\bar{\mu}>0$ need not be known, and we can build
proper models for which this assumption is satisfied.

######  Definition 2 (Model criticality.)

Let $\mathcal{M}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ be a model
satisfying Assumption 1. A point $\bar{x}\in X$ is said to be _model critical_
(M-critical for short) to problem (1) if

$\bar{x}\in\displaystyle\arg\min_{x\in
X}\mathcal{M}_{a}(x;\bar{x})\quad\forall\;a\in
A_{H}(\bar{x}):=\left\\{\alpha\in
A(\bar{x}):\,\mathcal{M}_{\alpha}(\bar{x};\bar{x})=H(\bar{x};\bar{x})\right\\}.$

Furthermore, $\bar{x}$ is said to be _feasible-model critical_ (FM-critical)
if it is M-critical and $c(\bar{x})\leq 0$. $\blacktriangleleft$

This new concept of criticality is a necessary optimality condition for
$\bar{x}\in X$ to solve problem (1).

###### Theorem 3.1

Let $H$ be the improvement function given in (2) and $\mathcal{M}$ be a model
satisfying Assumption 1. If $\bar{x}$ is a local solution to problem (1), then
$\bar{x}$ is FM-critical.

###### Proof

If $\bar{x}$ is a local solution to problem (1), then Lemma 1 (ii) asserts
that $\bar{x}\in\arg\min_{x\in X}H(x;\bar{x})$. Let $\bar{\mu}>0$ be as in
Assumption 1(b). Then,

$\displaystyle H(\bar{x};\bar{x})=\mathcal{M}(\bar{x};\bar{x})$
$\displaystyle\geq\displaystyle\min_{x\in
X}\;\mathcal{M}(x;\bar{x})+\frac{\bar{\mu}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}\geq\displaystyle\min_{x\in
X}\;H(x;\bar{x})=H(\bar{x};\bar{x}).$

Thus,

$\displaystyle H(\bar{x};\bar{x})$ $\displaystyle=\displaystyle\min_{x\in
X}\;\mathcal{M}(x;\bar{x})+\frac{\bar{\mu}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle=\displaystyle\min_{x\in X}\;\min_{a\in
A(\bar{x})}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\bar{\mu}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle\leq\displaystyle\min_{x\in X}\;\min_{a\in
A_{H}(\bar{x})}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\bar{\mu}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle=\displaystyle\min_{a\in A_{H}(\bar{x})}\;\min_{x\in
X}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\bar{\mu}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle\leq\displaystyle\min_{a\in
A_{H}(\bar{x})}\;\mathcal{M}_{a}(\bar{x};\bar{x})=H(\bar{x};\bar{x}),$

showing that $\bar{x}\in\arg\min_{x\in
X}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\bar{\mu}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
for all $a\in A_{H}(\bar{x})$. Convexity of $\mathcal{M}_{a}(\cdot;\bar{x})$
then implies that $\bar{x}\in\arg\min_{x\in X}\;\mathcal{M}_{a}(x;\bar{x})$
for all $a\in A_{H}(\bar{x})$, i.e., $\bar{x}$ is FM-critical. $\square$

##  4 Algorithm and convergence analysis

This section provides a proximal-like algorithm to compute M-critical points
under certain conditions. The approach’s working horse is a model
$\mathcal{M}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}$ satisfying
Assumption 1. We leave open the specific choice for the convex models
$\mathcal{M}_{a}$ composing $\mathcal{M}$, as we intend to make Algorithm 1
below applicable to a broad class of optimization problems of the form (1).
Under Assumption 1, we show that subproblem (3) below can be easily solved up
to a certain accuracy $\epsilon^{k}\geq 0$ at every iteration $k$. The
epigraphical nesting concept recalled in Section 2 then will be handy to show
that Algorithm 1 asymptotically computes an M-critical point. If the algorithm
starts with a feasible point, then it is ensured to compute an FM-critical
point.

Algorithm 1 Proximal-like Method

1:Let $x^{0}\in X$, $\kappa\in(0,1)$, $\mu^{0}\geq\kappa$, and
$\lambda\in[0,\kappa)$ be given

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

3: Let $y^{k}$ be an $\epsilon^{k}$-solution to  $\min_{y\in
X}\;{\mathcal{M}(y;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y-x^{k}\textstyle{\|}^{2},}\quad\qquad\mbox{with}\qquad\quad
0\leq\epsilon^{k}\leq\frac{\lambda}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
(3)

4: if $\textstyle{\|}y^{k}-x^{k}\textstyle{\|}\leq\texttt{Tol}$ then

5: Stop: and return $x^{k}$

6: end if

7: if $H(y^{k};x^{k})\leq
H(x^{k};x^{k})-\frac{\kappa-\lambda}{2}{\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}}$
then

8: Serious step: set $x^{k+1}=y^{k}$ and $\mu^{k+1}=\mu^{k}$

9: else

10: Null step: set $x^{k+1}=x^{k}$ and $\mu^{k+1}>\mu^{k}$

11: end if

12:end for

In Algorithm 1, the model is updated whenever a new serious step is performed:
in this case, the so-called stability center $x^{k}$ is changed in
$\mathcal{M}(\cdot;x^{k})$. For this task, the algorithm will need oracles
providing (partial) information on functions $f$ and $c$ (see equation (7) for
an example of how $x^{k}$ may impact $\mathcal{M}$).

Only $\mu^{k}$ changes after a null step. Assumption 1(b) is crucial to ensure
that once $\mu^{k}$ reaches a certain threshold $\bar{\mu}>0$, only serious
steps are performed. This is why we do not allow $\mu^{k}$ to decrease.

Line 3 requires solving inexactly a nonsmooth optimization problem, which we
assume has a solution (e.g. if $X$ is compact). The precision $\epsilon^{k}$
depends on the iterate $y^{k}$ and, thus, the optimization solver employed to
solve (3) does not receive $\epsilon^{k}$ as an input. However, this will not
be an issue in practice, as we explain in the following remark.

###### Remark 1

At every iteration $k$ subproblem (3) can be equivalently decomposed in a
finite family of convex optimization problems that can be solved by bundle
methods Correa_Lemarechal_1993 . As a result, we can globally solve the
nonconvex subproblem (3) by solving the following $|A(x^{k})|$ convex
subproblems

$y_{a}^{k}:=\displaystyle\arg\min_{x\in
X}\mathcal{M}_{a}(x;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2},$
(4)

and letting $y^{k}:=y^{k}_{a^{*}}$, with $a^{*}$ the index yielding the
minimum functional value above. In Appendix A we rely on
Correa_Lemarechal_1993  to detail a bundle-like algorithm to solve the above
subproblems up to tolerance $\epsilon^{k}\geq 0$ satisfying the condition on
the right of equation (3).

###### Theorem 4.1

Consider problem (1), a model $\mathcal{M}$ satisfying Assumption 1, and
suppose that subproblem (3) is solvable for all $k$. Let $\\{x^{k}\\}_{k}$ be
the sequence produced by Algorithm 1 with $\texttt{Tol}=0$. Furthermore,
suppose that

  1. (i)

$\\{x^{k}\\}_{k}$ is bounded;

  2. (ii)

after a null step, the rule to update the prox-parameter ensures that
$\mu^{k+1}>\mu^{k}+\delta$ for some constant $\delta>0$;

  3. (iii)

for any cluster point $\bar{x}$ of $\\{x^{k}\\}_{k}$, there exists a
subsequence $\\{x^{k_{i}}\\}_{i}$ such that $\lim_{i}x^{k_{i}}=\bar{x}$ and
$\operatorname{epi}\left(\bar{\mathcal{M}}(\cdot;\bar{x})+{\rm
i}_{X}\right)\subseteq\operatorname{LimInn}\left(\operatorname{epi}\left(\mathcal{M}(\cdot;x^{k_{i}})+{\rm
i}_{X}\right)\right)$, with $\bar{\mathcal{M}}(\cdot;\bar{x}):=\min_{a\in
A_{H}(\bar{x})}\mathcal{M}_{a}(\cdot;\bar{x})$.

Then $\bar{x}$ is M-critical, which upgrades to FM-critical if $\bar{x}$ is
feasible (that is certainly true if the algorithm starts with a feasible
point).

###### Proof

First, observe that if the algorithm stops at iteration $k$ with
$\texttt{Tol}=0$, then $y^{k}=x^{k}$ and, thus, condition (3) implies that
$\bar{x}=x^{k}$ is an exact solution to subproblem (3). Thus,

$\displaystyle H(\bar{x};\bar{x})$
$\displaystyle=\mathcal{M}(\bar{x};\bar{x})=\min_{x\in
X}\;\mathcal{M}(x;\bar{x})+\frac{\mu^{k}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle\leq\min_{x\in X}\;\min_{a\in
A_{H}(\bar{x})}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\mu^{k}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle\leq\min_{x\in
X}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\mu^{k}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}\quad\forall\;a\in
A_{H}(\bar{x})$ $\displaystyle\leq H(\bar{x};\bar{x}),$

showing that $\bar{x}\in\arg\min_{x\in
X}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\mu^{k}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
for all $a\in A_{H}(\bar{x})$. Convexity of $\mathcal{M}_{a}(\cdot;x)$ yields
$\bar{x}\in\arg\min_{x\in X}\;\mathcal{M}_{a}(x;\bar{x})$ for all $a\in
A_{H}(\bar{x})$, i.e., $\bar{x}$ is M-critical.

From now on, let us assume that the algorithm loops forever:

Claim 1: Algorithm 1 produces finitely many null steps. Note that if
$\mu^{k}\geq\bar{\mu}+\kappa$ (with $\bar{\mu}$ ensured by Assumption 1), then
$y^{k}$ is declared a serious iterate. Indeed,

$\displaystyle H(y^{k};x^{k})$
$\displaystyle\leq\mathcal{M}(y^{k};x^{k})+\frac{\bar{\mu}}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
$\displaystyle=\mathcal{M}(y^{k};x^{k})+\frac{\bar{\mu}+\kappa}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}-\frac{\kappa}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
$\displaystyle\leq\min_{x\in
X}\left\\{\mathcal{M}(x;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}\right\\}+\epsilon^{k}-\frac{\kappa}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
$\displaystyle\leq\mathcal{M}(x^{k};x^{k})+\epsilon^{k}-\frac{\kappa}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
$\displaystyle\leq
H(x^{k};x^{k})-\frac{\kappa-\lambda}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2},$

where the first inequality is due to by Assumption 1(b), and the two last ones
follow from definition of $y^{k}$ and identity
$\mathcal{M}(x^{k};x^{k})=H(x^{k};x^{k})$. Thus $y^{k}$ satisfies the descent
test and becomes a serious iterate. As $\mu^{k}$ increases by at least
$\delta>0$ after a null step (hypothesis (ii)), we conclude that Algorithm 1
generates only finitely many null steps. Hence, there exists $\bar{k}$ such
that, for all $k\geq\bar{k}$,

$\mu^{k}=\mu^{\prime}:=\mu^{\bar{k}},\quad
x^{k+1}\in\epsilon^{k}\mbox{-}\operatorname{argmin}_{x\in
X}\mathcal{M}(x;x^{k})+\frac{\mu^{\prime}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2},\quad\mbox{and}$
$H(x^{k+1};x^{k})\leq
H(x^{k};x^{k})-\frac{\kappa-\lambda}{2}\textstyle{\|}x^{k+1}-x^{k}\textstyle{\|}^{2}.$

Claim 2: The sequence $\\{\textstyle{\|}y^{k}-x^{k}\textstyle{\|}\\}_{k}$
vanishes. First, note that $x^{k+1}=y^{k}$ for all $k\geq\bar{k}$. We split
our analysis into two cases: there exists an iteration counter
$k_{1}\geq\bar{k}$ such that $x^{k_{1}}$ is feasible for problem (1) or such a
$k_{1}$ does not exist. If we assume the existence of such $k_{1}$, by an
inductive argument on Lemma 1 (item (i.1) with $\mu=\kappa-\lambda>0$) we can
deduce that for every $k\geq k_{1}$ the iterate $x^{k}\in X^{c}$ and the
following inequality holds

$f(x^{k+1})\leq
f(x^{k})-\frac{\kappa-\lambda}{2}\textstyle{\|}x^{k+1}-x^{k}\textstyle{\|}^{2}.$

Rearranging the terms and adding the last expression for every $k\geq k_{1}$,
we obtain

$\sum_{k=k_{1}}^{\infty}\textstyle{\|}x^{k+1}-x^{k}\textstyle{\|}^{2}\leq\frac{2}{\kappa-\lambda}\sum_{k=k_{1}}^{\infty}\left(f(x^{k+1})-f(x^{k})\right)\leq\frac{2}{\kappa-\lambda}\left(f(x^{k_{1}})-\lim_{k\to\infty}f(x^{k}))\right).$

It follows by the continuity of $f$, and boundedness of $\\{x^{k}\\}_{k}$ that
the right-hand side of the above inequality is finite, and so is the sum at
the left-hand side. This implies that
$\\{\textstyle{\|}x^{k+1}-x^{k}\textstyle{\|}\\}_{k}$ converges to zero. If we
now assume that $x^{k}$ is not feasible for every $k\geq\bar{k}$ then,
inequality

$c(x^{k+1})\leq
c(x^{k})-\frac{\kappa-\lambda}{2}\textstyle{\|}x^{k+1}-x^{k}\textstyle{\|}^{2},$

holds for every $k\geq\bar{k}$ by applying induction on Lemma 1 (item (i.2)).
Then, by a similar argument as the one used in the previous case, we deduce
that the sequence $\\{\textstyle{\|}x^{k+1}-x^{k}\textstyle{\|}\\}_{k}$
vanishes. Hence, $\epsilon^{k}\to 0$ (due to (3)) and $\\{y^{k_{i}}\\}_{i}$
also converges to $\bar{x}$:

$\lim_{i}y^{k_{i}}=\lim_{i}x^{k_{i}+1}=\lim_{i}x^{k_{i}}=\bar{x}.$

Claim 3: Epigraphical nesting of the sequence. Remark that, using hypothesis
(iii), we can consider a sequence $z^{i}\to x$ for any given
$x\in\mathbb{R}^{n}$ holding item (ii) of Lemma 2 for
$\mathcal{M}(\cdot;x^{k_{i}})+{\rm i}_{X}$. Thus,

$\limsup_{i}\mathcal{M}(z^{i};x^{k_{i}})+{\rm
i}_{X}(z^{i})+\frac{\mu^{k_{i}}}{2}\textstyle{\|}z^{i}-x^{k_{i}}\textstyle{\|}^{2}\leq\bar{\mathcal{M}}(x;\bar{x})+{\rm
i}_{X}(x)+\frac{\mu^{\prime}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$

which, again by item (ii) in Lemma 2, leads to

$\operatorname{epi}\left(\bar{\mathcal{M}}(\cdot;\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}\cdot-\bar{x}\textstyle{\|}^{2}+{\rm
i}_{X}\right)\subseteq\operatorname{LimInn}\left(\operatorname{epi}\left(\mathcal{M}(\cdot;x^{k_{i}})+\frac{\mu^{k_{i}}}{2}\textstyle{\|}\cdot-x^{k_{i}}\textstyle{\|}^{2}+{\rm
i}_{X}\right)\right).$

Observe that, for all $k$,

$\displaystyle H(y^{k};x^{k})$
$\displaystyle\leq\mathcal{M}(y^{k};x^{k})+\frac{\bar{\mu}}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
$\displaystyle=\mathcal{M}(y^{k};x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}+\frac{\bar{\mu}-\mu^{k}}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
$\displaystyle\leq\min_{y\in
X}\mathcal{M}(y;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y-x^{k}\textstyle{\|}^{2}+\epsilon^{k}+\frac{\bar{\mu}-\mu^{k}}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}$
$\displaystyle\leq
H(x^{k};x^{k})+\frac{\bar{\mu}-\mu^{k}+\lambda}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}.$

Thus, for the given subsequences and taking the limit with $i$ going to
infinity, the continuity of $H$ leads to

$\displaystyle\lim_{i}\mathcal{M}(y^{k_{i}};x^{k_{i}})+\frac{\mu^{k}}{2}\textstyle{\|}y^{k_{i}}-x^{k_{i}}\textstyle{\|}^{2}+{\rm
i}_{X}(y^{k_{i}})$ $\displaystyle=H(\bar{x};\bar{x})$
$\displaystyle=\mathcal{M}_{a}(\bar{x};\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}\bar{x}-\bar{x}\textstyle{\|}^{2}\quad\forall\;a\in
A_{H}(\bar{x})$
$\displaystyle=\bar{\mathcal{M}}(\bar{x};\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}\bar{x}-\bar{x}\textstyle{\|}^{2}.$

We can thus conclude that $\bar{x}\in\operatorname{argmin}_{x\in
X}\bar{\mathcal{M}}(x;\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
thanks to Claims 2-3 and (RoysetWets_2022, , Prop. 6.14), which in turn gives

$\displaystyle H(\bar{x};\bar{x})$ $\displaystyle=\min_{x\in
X}\bar{\mathcal{M}}(x;\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle=\min_{x\in X}\;\min_{a\in
A_{H}(\bar{x})}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
$\displaystyle\leq\min_{x\in
X}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}\quad\forall\;a\in
A_{H}(\bar{x})$ $\displaystyle\leq H(\bar{x};\bar{x}),$

showing that $\bar{x}\in\arg\min_{x\in
X}\;\mathcal{M}_{a}(x;\bar{x})+\frac{\mu^{\prime}}{2}\textstyle{\|}x-\bar{x}\textstyle{\|}^{2}$
for all $a\in A_{H}(\bar{x})$. Convexity of $\mathcal{M}_{a}(\cdot;x)$ yields
$\bar{x}\in\arg\min_{x\in X}\;\mathcal{M}_{a}(x;\bar{x})$ for all $a\in
A_{H}(\bar{x})$, i.e., $\bar{x}$ is M-critical. In addition, if $x^{0}\in
X^{c}=\left\\{x\in X:\,c(x)\leq 0\right\\}$, then the algorithm ensures that
$\\{x^{k}\\}_{k}\subset X^{c}$ due to Lemma 1 (i.1). FM-criticality of
$\bar{x}$ follows from the fact that $X^{c}$ is closed (as $X$ is closed and
$c$ is continuous). $\square$

We conclude this section with a few comments on the hypothesis employed by
Theorem 4.1. Note that hypothesis (i) is simple to satisfy (e.g. if $X$ is
compact). Hypothesis (ii) is a trivial one. Its goal is to ensure that if the
algorithm produces infinitely many null steps, then $\mu^{k}\to\infty$. The
more stringent conditions are the ones related to the model, i.e., Assumption
1 and hypothesis (iii). We recall that the threshold $\bar{\mu}>0$ in
Assumption 1 is not required to be known. In the next sections, we show how to
build proper models for which Assumption 1 and hypothesis (iii) are satisfied.

##  5 Optimization problems from stochastic programming

The following structure is motivated by reliability-based optimization
problems, where the design cost of an engineering system is minimized while
satisfying reliability constraints Byun_Oliveira_Royset_2023 . The structure
covers chance constraints and buffered optimization models constructed via
sample average approximation. Details are postponed to the section on
numerical experiments. Let us focus on the resulting deterministic
optimization problem, which is a particular case of problem (1)

$\left\\{\begin{array}[]{lll}\displaystyle\min_{x\in
X}&f(x):=\displaystyle\sum_{j=1}^{N_{f}}\max_{\ell\in
F_{j}}\\{\mathfrak{f}^{1}_{j\ell}(x)+\mathfrak{f}^{2}_{j\ell}(x)\\}\\\
\mbox{s.t.}&c(x):=\displaystyle\sum_{j=1}^{N_{c}}\max_{\ell\in
C_{j}}\\{\mathfrak{c}^{1}_{j\ell}(x)+\mathfrak{c}^{2}_{j\ell}(x)\\}\leq
0,\end{array}\right.$ (5)

under the following assumptions:

  1. (a)

$\emptyset\neq X$ is a convex and closed subset of $\mathbb{R}^{n}$;
$\mathfrak{f}^{1}_{j\ell},\mathfrak{c}^{1}_{j\ell}:\mathbb{R}^{n}\to\mathbb{R}$
are convex functions,
$\mathfrak{f}^{2}_{j\ell},\mathfrak{c}^{2}_{j\ell}:\mathbb{R}^{n}\to\mathbb{R}$
are weakly-concave functions, and $F_{j}$, $C_{j}$ are finite index sets.

Recall that a function $\phi:\mathbb{R}\to\mathbb{R}$ is weakly-concave if
there exists a modulus $\mu_{\phi}\geq 0$ such that
$\phi(x)-\frac{\mu}{2}\textstyle{\|}x\textstyle{\|}^{2}$ is concave for all
$\mu\geq\mu_{\phi}$. This implies that, for all $\mu\geq\mu_{\phi}$ and all
$\xi_{\phi}^{x}\in\partial^{\tt C}\phi(x)$, the following inequality holds

$\phi(y)\leq\phi(x)+\langle\xi_{\phi}^{x},y-x\rangle+\frac{\mu}{2}\textstyle{\|}y-x\textstyle{\|}^{2}\quad\forall\,y\in\mathbb{R}^{n}.$
(6)

In what follows we work with the improvement function $H$ of this problem and
a model $\mathcal{M}$ for $H$ satisfying Assumption 1 and hypothesis (iii) in
Theorem 4.1. We will work with Clarke subdifferentials, which are well defined
since the finite-valued functions $f$ and $c$ in (5) are locally Lipschitz.
Our approach is general in the sense that we do not assume available oracles
for $f$ and $c$. Instead, we assume that there are individual oracles for the
component functions in (5), providing us with function values and at least one
subgradient for every component function: given $x\in X$, the oracles provide
us with $\mathfrak{f}^{2}_{j\ell}(x)$, $\zeta_{{j\ell}}\in\partial^{\tt
C}\mathfrak{f}^{2}_{j\ell}(x)$, $\mathfrak{c}^{2}_{j\ell}(x)$, and
$\xi_{{j\ell}}\in\partial^{\tt C}\mathfrak{c}^{2}_{j\ell}(x)$. With this
partial information, let us define the following convex model:

$\mathcal{M}_{a}(y;x):=\max\left\\{\begin{array}[]{ll}\displaystyle\sum_{j=1}^{N_{f}}\max_{\ell\in
F_{j}}\left\\{\mathfrak{f}^{1}_{j\ell}(y)+\mathfrak{f}^{2}_{j\ell}(x)+\langle\zeta_{{j\ell}},y-x\rangle\right\\}-\tau_{f}(x),\\\
\displaystyle\sum_{j=1}^{N_{c}}\max_{\ell\in
C_{j}}\left\\{\mathfrak{c}^{1}_{j\ell}(y)+\mathfrak{c}^{2}_{j\ell}(x)+\langle\xi_{{j\ell}},y-x\rangle\right\\}\end{array}\right\\}$
(7a) where $a$ denotes the tuple $(\zeta,\xi)$ with
$\zeta=(\zeta_{j\ell})_{j\ell}$ and $\xi=(\xi_{j\ell})_{j\ell}$ being the
subgradients given by the oracles. If the oracles can provide us with more
than one tuple, say a set
$A(x)\subseteq\left\\{(\zeta,\xi):\zeta_{{j\ell}}\in\partial^{\tt
C}\mathfrak{f}^{2}_{j\ell}(x)\;,\;\xi_{{j\ell}}\in\partial^{\tt
C}\mathfrak{c}^{2}_{j\ell}(x)\quad\forall\,j,\ell\right\\}$ (7b) of them, we
may consider the model $\mathcal{M}(y;x)=\min_{a\in
A(x)}\;\mathcal{M}_{a}(y;x).$ (7c)

Observe that $\mathcal{M}_{a}(\cdot;x)$ is convex and
$\mathcal{M}_{a}(x;x)=H(x;x)$ for all $x$ and $a\in A(x)$, yielding thus
Assumption 1(a). In the remainder of this section we assume that

$|A(x)|\leq{\tt n}<\infty,\quad\mbox{with ${\tt n}\geq 1$ an integer number.}$
(8)

This means that we get at most ${\tt n}$ subgradient tuples from the oracles.
In what follows we proceed to show that $\mathcal{M}$ above satisfies also
Assumption 1(b) and hypothesis (iii) in Theorem 4.1.

###  5.1 Checking the required assumptions on the model

The following result shows that model $\mathcal{M}$ in (7) satisfies
Assumption 1(b).

###### Lemma 3

Let each function
$\mathfrak{f}^{2}_{j\ell},\mathfrak{c}^{2}_{j\ell}:\mathbb{R}^{n}\to\mathbb{R}$
in (5) be weakly concave. Then there exists $\bar{\mu}\geq 0$ such, for all
$\mu\geq\bar{\mu}$,

$H(y;x)\leq\mathcal{M}(y;x)+\frac{\mu}{2}\textstyle{\|}y-x\textstyle{\|}^{2}\quad\forall\;y,x\in\mathbb{R}^{n}.$

###### Proof

Let $\mu_{j\ell}\geq 0$ be the modulus of the weakly-concave function
$\mathfrak{f}^{2}_{j\ell}$. If follows from (6) that
$\mathfrak{f}^{2}_{j\ell}(y)\leq\mathfrak{f}^{2}_{j\ell}+\langle\zeta_{{j\ell}},y-x\rangle+\frac{\mu}{2}\textstyle{\|}y-x\textstyle{\|}^{2}$
for all $\mu\geq\mu_{j\ell}$, independently of the choice of $\zeta$. As every
index set $F_{i}$ in (5) are finite and $j=1,\ldots,N_{f}$, the value
$\bar{\mu}_{f}:=\max_{j=1,\ldots,N_{f}}\max_{\ell\in F_{j}}\mu_{j\ell}$ is
finite. Hence, for all $y\in\mathbb{R}^{n}$ and all $\mu\geq\bar{\mu}_{f}$,

$\displaystyle f(y)$ $\displaystyle\leq\sum_{j=1}^{N_{f}}\max_{\ell\in
F_{j}}\left\\{\mathfrak{f}^{1}_{j\ell}(y)+\mathfrak{f}^{2}_{j\ell}(x)+\langle\zeta_{{j\ell}},y-x\rangle+\frac{\bar{\mu}_{f}}{2}\textstyle{\|}y-x\textstyle{\|}^{2}\right\\}$
$\displaystyle=\sum_{j=1}^{N_{f}}\max_{\ell\in
F_{j}}\left\\{\mathfrak{f}^{1}_{j\ell}(y)+\mathfrak{f}^{2}_{j\ell}(x)+\langle\zeta_{{j\ell}},y-x\rangle\right\\}+\frac{N_{f}\bar{\mu}_{f}}{2}\textstyle{\|}y-x\textstyle{\|}^{2}.$

An analogous argument concerning functions $\mathfrak{c}^{2}_{j\ell}$ ensures
the existence of a threshold $\bar{\mu}_{c}<\infty$ such that

$c(y)\leq\sum_{j=1}^{N_{c}}\max_{\ell\in
C_{j}}\left\\{\mathfrak{c}^{1}_{j\ell}(y)+\mathfrak{c}^{2}_{j\ell}(x)+\langle\xi_{{j\ell}},y-x\rangle\right\\}+\frac{N_{c}\bar{\mu}_{c}}{2}\textstyle{\|}y-x\textstyle{\|}^{2}.$

Then, by taking
$\bar{\mu}=\max\left\\{N_{f}\bar{\mu}_{f},N_{c}\bar{\mu}_{c}\right\\}$ and
recalling the definitions of $H$ in (2) and $\mathcal{M}_{a}$ in (7a), we
conclude that
$H(y;x)\leq\mathcal{M}_{a}(y;x)+\frac{\bar{\mu}}{2}\textstyle{\|}y-x\textstyle{\|}^{2}$.
The result follows by taking the minimum over $a\in A(x)$ given in (7b) and
recalling (7c). $\square$

Next, we show that $\mathcal{M}$ in (7) also satisfies hypothesis (iii) in
Theorem 4.1. To that end, we shall keep in mind the following observation.

###### Remark 2

For any given cluster point $\bar{x}$, there can be several subsequences of
subgradients of the same function $\mathfrak{f}^{2}_{j\ell}$ converging to
different cluster points, all belonging to $\partial^{\tt
C}\mathfrak{f}^{2}_{j\ell}(\bar{x})$ (analogously for the function
$\mathfrak{c}^{2}_{j\ell}$). As a result, there can be several different
models $\bar{\mathcal{M}}(\cdot;\bar{x})$, any of which suffices to show that
hypothesis (iii) in Theorem 4.1 is satisfied. Indeed, recall that given a
bounded sequence $\\{x^{k}\\}_{k}$, then the sequences
$\\{\zeta_{{j\ell}}^{k}\\}_{k}$ and $\\{\xi_{{j\ell}}^{k}\\}_{k}$ of
subgradients are bounded as well. Given an arbitrary cluster point $\bar{x}$
of $\\{x^{k}\\}_{k}$, we can therefore extract convergent subsequences,
$\\{x^{k_{i}}\\}_{i}$, $\\{\zeta_{{j\ell}}^{k_{i}}\\}_{i}$ and
$\\{\xi_{{j\ell}}^{k_{i}}\\}_{i}$ such that

$\lim_{i}x^{k_{i}}=\bar{x},\quad\lim_{i}\zeta_{{j\ell}}^{k_{i}}=\bar{\zeta}_{{j\ell}}\in\partial^{\tt
C}\mathfrak{f}^{2}_{j\ell}(\bar{x}),\quad\mbox{and}\quad\lim_{i}\xi_{{j\ell}}^{k_{i}}=\bar{\xi}_{{j\ell}}\in\partial^{\tt
C}\mathfrak{c}^{2}_{j\ell}(\bar{x}),\quad\forall j,\ell$ (9)

because the Clarke subdifferential is locally bounded and outer-semicontinuous
(Clarke_1987, , Prop. 2.1.5). This limiting tuple of subgradients defines
$\bar{a}$: under assumption (8), such tuple is associated with some
$a^{k_{i}^{\prime}}\in A(x^{k_{i}^{\prime}})$ (in the sense that
$\lim_{i}a^{k_{i}^{\prime}}\to\bar{a}$) and thus we can define $A(\bar{x})$ as
$\operatorname{LimOut}\left(A(x^{k_{i}})\right)$. Furthermore, as each index
$a\in A(x^{k})$ satisfies $\mathcal{M}_{a}(x^{k};x^{k})=H(x^{k};x^{k})$, we
have that $A_{H}(x^{k})=A(x^{k})$ for all $k$. All in all, we can conclude
that

$A_{H}(\bar{x})=A(\bar{x}):=\operatorname{LimOut}\left(A(x^{k_{i}})\right).$

$\blacktriangleleft$

Hence, we can define the following model at cluster point $\bar{x}$ by
employing those cluster points of subgradients defined in (9):

$\displaystyle\bar{\mathcal{M}}(y;\bar{x})$ $\displaystyle:=\min_{a\in
A_{H}(\bar{x})}\mathcal{M}_{a}(y;\bar{x})$ $\displaystyle=\min_{a\in
A_{H}(\bar{x})}\max\left\\{\begin{array}[]{ll}\displaystyle\sum_{j=1}^{N_{f}}\max_{\ell\in
F_{j}}\left\\{\mathfrak{f}^{1}_{j\ell}(y)+\mathfrak{f}^{2}_{j\ell}(\bar{x})+\langle\bar{\zeta}_{{j\ell}},y-\bar{x}\rangle\right\\}-\tau_{f}(\bar{x}),\\\
\displaystyle\sum_{j=1}^{N_{c}}\max_{\ell\in
C_{j}}\left\\{\mathfrak{c}^{1}_{j\ell}(y)+\mathfrak{c}^{2}_{j\ell}(\bar{x})+\langle\bar{\xi}_{{j\ell}},y-\bar{x}\rangle\right\\}\end{array}\right\\}.$
(12)

###### Theorem 5.1

Given model $\mathcal{M}$ in (7) satisfying (8), assume hypotheses (i) and
(ii) in Theorem 4.1 hold. Then every cluster point $\bar{x}$ of the sequence
$\\{x^{k}\\}_{k}$ computed by Algorithm 1 is M-critical. Furthermore, provided
that an iterate $x^{k}$ is feasible, $\bar{x}$ is FM-critical.

Proof Given the subsequence $\\{x^{k_{i}}\\}_{i}$ from (9) and
$\bar{\mathcal{M}}(\cdot;\bar{x})$ as in (12), the result follows from Theorem
4.1 once show that
$\operatorname{epi}\left(\bar{\mathcal{M}}(\cdot;\bar{x})+{\rm
i}_{X}\right)\subseteq\operatorname{LimInn}\left(\operatorname{epi}\left(\mathcal{M}(\cdot;x^{k_{i}})+{\rm
i}_{X}\right)\right)$. To this end, consider an arbitrary point of $y\in X$
and any sequence $\\{y^{i}\\}_{i}\subseteq X$ such that $y^{i}\to y$. For the
specific sequences in (9) we get (see Remark 2)

$\displaystyle\bar{\mathcal{M}}(y;\bar{x})$ $\displaystyle=\min_{a\in
A_{H}(\bar{x})}\mathcal{M}_{a}(y;\bar{x})=\min_{a\in\operatorname{LimOut}\left(A(x^{k_{i}})\right)}\mathcal{M}_{a}(y;\bar{x})$
$\displaystyle=\lim_{i}\min_{a\in
A(x^{k_{i}})}\mathcal{M}_{a}(y^{i};x^{k_{i}})=\lim_{i}\mathcal{M}(y^{i};x^{k_{i}}),$

where the third inequality follows from our construction of
$\operatorname{LimOut}\left(A(x^{k_{i}})\right)$ and continuity of
$\mathcal{M}_{a}(\cdot;\cdot)$ on both arguments. Lemma 2 item (ii) yields the
wished inclusion. $\square$

We stress that in the above theorem, the only theoretical gain in working with
$A(x^{k})>1$ (i.e., a nonconvex model) is the fact that $A_{H}(\bar{x})$ in
our criticality definition 2 can have more than one index. Thus, the
criticality we get in practice might be more robust: a richer model will
likely yield better results. In Section 6, we consider a specific DC setting
and show that a richer model yields $B$-stationary, the strongest stationary
condition in DC programming.

###  5.2 Particular instances

Let us now consider simpler variants of problem (5).

#### Case 1:

Consider problem (5) with $N_{f}=N_{c}=1$ and $F_{1}=C_{1}=\\{1\\}$. In this
case, $f(x)=f_{1}(x)+f_{2}(x)$, $f(x)=c_{1}(x)+c_{2}(x)$ with $f_{1},c_{1}$
convex and $f_{2},c_{2}$ weakly-concave functions. Moreover, let us consider
$f_{2}$ and $c_{2}$ to be smooth (i.e., continuously differentiable). Model
(7) reduces to

$\mathcal{M}(y;x)=\max\left\\{\begin{array}[]{l}f_{1}(y)+f_{2}(x)+\langle\nabla
f_{2}(x),y-x\rangle-\tau_{f}(x)\\\ c_{1}(y)+c_{2}(x)+\langle\nabla
c_{2}(x),y-x\rangle\end{array}\right\\}$

as $A(x)=\\{(\nabla f_{2}(x),\nabla c_{2}(x))\\}$. If $\bar{x}\in X$ is FM-
critical for the model
$\bar{\mathcal{M}}(\cdot;\bar{x})=\mathcal{M}(\cdot;\bar{x})$ (which coincides
with the one in (12) due to continuity of the gradients), we would have that

$\bar{x}\in\arg\min_{x\in X}\;\mathcal{M}(x;\bar{x});\quad 0\geq
c(\bar{x})=[c_{1}(\cdot)+c_{2}(\bar{x})+\langle\nabla
c_{2}(x),\cdot-\bar{x}\rangle](\bar{x}).$

In that case, and assuming the Constraint Qualification (CQ)
$c(\bar{x})=0\Longrightarrow 0\not\in\partial c_{1}(\bar{x})+\nabla
c_{2}(\bar{x})+N_{X}(\bar{x})$ holds, there exists $\bar{u}\in\mathbb{R}$ such
that the following generalized KKT system is satisfied

$\left\\{\begin{array}[]{llllll}0\in\partial f_{1}(\bar{x})+\nabla
f_{2}(\bar{x})+\bar{u}[\partial c_{1}(\bar{x})+\nabla
c_{2}(\bar{x})]+N_{X}(\bar{x})\\\ \bar{x}\in X,\;c(\bar{x})\leq
0,\;\bar{u}\geq 0,\;\bar{u}c(\bar{x})=0.\end{array}\right.$

Note that, as $\partial f_{1}(\bar{x})+\nabla f_{2}(\bar{x})=\partial^{\tt
C}f(\bar{x})$ and $\partial c_{1}(\bar{x})+\nabla c_{2}(\bar{x})=\partial^{\tt
C}c(\bar{x})$ (Clarke_1987, , Prop. 2.3.3), our FM-criticality boils down to
the standard KKT optimality conditions (stationarity) in nonsmooth nonlinear
programming, given the stated CQ holds.

In the setting of convex plus weakly-concave functions a bundle method with
convergence guarantees was recently proposed in Ksenia_2023 .

#### Case 2 (DC setting):

Here, not only $f_{1}$ and $c_{1}$ but also $-f_{2}$ and $-c_{2}$ are assumed
to be convex functions. Accordingly, the threshold $\bar{\mu}$ in Assumption 1
is zero. Furthermore, the criticality condition is as in the previous
paragraph, replacing $\nabla f_{2}$ and $\nabla c_{2}$ by $\partial^{\tt
C}f_{2}$ and $\partial^{\tt C}c_{2}$, boiling down to the usual criticality
condition in DC-constrained DC programming Pang_2017 ; Javal_2021 . As we
discuss next, stronger conditions are possible after adding assumptions and
making changes in the model $\mathcal{M}$.

##  6 Special DC setting

We have just shown that DC-constrained DC programming problems fit into
formulation (5), and thus the M-criticality coincides (when ${\tt n}=1$ in
(8)) with the usual criticality concept employed in the DC literature. The
goal here is to consider a special class of DC problems, a different model,
and show that Algorithm 1 obtains the sharpest stationarity definition in DC
programming, i.e., $B$-stationarity. The family of DC programs we investigate
is the one considered in Pang_2017 : problem (1) with

$f(x):=f_{1}(x)+f_{2}(x)$, $c(x)=c_{1}(x)+c_{2}(x)$,
$f_{1},c_{1}:\mathbb{R}^{n}\to\mathbb{R}$ convex, and  (13a)
$f_{2}(x):=\min_{j=1,\dots,m_{f}}\mathfrak{f}_{j}(x)\quad\text{and}\quad
c_{2}(x):=\min_{\ell=1,\dots,m_{c}}\mathfrak{c}_{\ell}(x),$ (13b)

where all $\mathfrak{f}_{j},\mathfrak{c}_{\ell}:\mathbb{R}^{n}\to\mathbb{R}$
are continuously differentiable concave functions. For $\varepsilon\geq 0$,
let us define the following index set:

$A_{f}^{\varepsilon}(x):=\left\\{j:\mathfrak{f}_{j}(x)\leq
f_{2}(x)+\varepsilon\right\\},\quad\mbox{and}\quad
A_{c}^{\varepsilon}(x):=\left\\{\ell:\mathfrak{c}_{\ell}(x)\leq
c_{2}(x)+\varepsilon\right\\}.$

For $\varepsilon=0$, these sets reduce to the active index sets of each
function. Let $\bar{x}\in X^{c}$ be a point satisfying the following
_pointwise Slater_ CQ: there exist $\ell\in A_{c}^{0}(\bar{x})$, and
$d^{\ell}\in\mathcal{T}_{X}(\bar{x})$ such that
$c^{\prime}_{1}(\bar{x};d^{\ell})<\langle\nabla\mathfrak{c}_{\ell}(\bar{x}),d^{\ell}\rangle$.
Under this CQ, it was shown in (Ackooij_Oliveira_2022, , Lem. 2.1) that
$\bar{x}$ is an $B$-stationary point to problem (5) under the setting of (13)
if and and only if $\bar{x}$ solves the following convex subproblems for all
$j\in A_{f}^{0}(\bar{x})$ and $\ell\in A_{c}^{0}(\bar{x})$ (see also
(Pang_2017, , Lem. 1) for a somehow stronger requirement):

$\left\\{\begin{array}[]{lll}\displaystyle\min_{x\in
X}&f_{1}(x)+\mathfrak{f}_{j}(\bar{x})+\langle\nabla\mathfrak{f}_{j}(\bar{x}),x-\bar{x}\rangle\\\
\mbox{s.t.}&c_{1}(x)+\mathfrak{c}_{\ell}(\bar{x})+\langle\nabla\mathfrak{c}_{\ell}(\bar{x}),x-\bar{x}\rangle\leq
0.\end{array}\right.$ (14)

With this in mind, we propose to furnish Algorithm 1 with the following model
for $H$:

$\mathcal{M}(y;x)=\max\left\\{\begin{array}[]{l}f_{1}(y)+\displaystyle\min_{j\in
A_{f}^{\varepsilon}(x)}[\mathfrak{f}_{j}(x)+\langle\nabla\mathfrak{f}_{j}(x),y-x\rangle]-\tau_{f}(x)\\\
c_{1}(y)+\displaystyle\min_{\ell\in
A_{c}^{\varepsilon}(x)}[\mathfrak{c}_{2}(x)+\langle\nabla\mathfrak{c}_{\ell}(x),y-x\rangle]\end{array}\right\\}.$
(15)

This nonconvex model can be equivalently written as

$\displaystyle\mathcal{M}(y;x)$ $\displaystyle=\displaystyle\min_{j\in
A_{f}^{\varepsilon}(x)}\displaystyle\min_{\ell\in
A_{c}^{\varepsilon}(x)}\max\left\\{\begin{array}[]{l}f_{1}(y)+[\mathfrak{f}_{j}(x)+\langle\nabla\mathfrak{f}_{j}(x),y-x\rangle]-\tau_{f}(x)\\\
c_{1}(y)+[\mathfrak{c}_{2}(x)+\langle\nabla\mathfrak{c}_{\ell}(x),y-x\rangle]\end{array}\right\\}$
(18) $\displaystyle=\displaystyle\min_{j\in
A_{f}^{\varepsilon}(x)}\displaystyle\min_{\ell\in
A_{c}^{\varepsilon}(x)}\mathcal{M}_{j\ell}(y;x),$ (19)

with every $\mathcal{M}_{j\ell}(\cdot;x)$ convex. Thus, we are in the setting
of Assumption 1 with $\bar{\mu}=0$ and $A(x)=A^{\varepsilon}_{f}(x)\times
A^{\varepsilon}_{c}(x)$.

###### Lemma 4

Let $\varepsilon\geq 0$ and suppose that $\bar{x}$ is FM-critical.
Furthermore, suppose that the pointwise Slater CQ is satisfied if
$c(\bar{x})=0$. Then $\bar{x}$ is an $B$-stationary point to problem (1)
having structure (13).

###### Proof

Note that model (19) satisfies, for all $\varepsilon\geq 0$, Assumption 1 (the
threshold $\bar{\mu}$ is zero in the DC setting). Recalling Remark 1 and
assuming $\bar{x}$ to be FM-critical

$\bar{x}\in\displaystyle\arg\min_{x\in
X}\;\mathcal{M}_{j\ell}(x;\bar{x})\quad\forall\,j\in
A_{f}^{0}(\bar{x}),\;\forall\,\ell\in A_{x}^{0}(\bar{x})$

and $c(\bar{x})\leq 0$, yielding
$\mathcal{M}_{j\ell}(\bar{x};\bar{x})=H(\bar{x};\bar{x})=0$. Hence,

$f_{1}(x)+\mathfrak{f}_{j}(\bar{x})+\langle\nabla\mathfrak{f}_{j}(\bar{x}),x-\bar{x}\rangle\geq
f(\bar{x})\quad\forall\;x\in
X:\;c_{1}(x)+\mathfrak{c}_{\ell}(\bar{x})+\langle\nabla\mathfrak{c}_{\ell}(\bar{x}),x-\bar{x}\rangle\leq
0$

and all $j\in A_{f}^{0}(\bar{x})$ and $\ell\in A_{x}^{0}(\bar{x})$. In
particular, $\bar{x}$ solves

$\left\\{\begin{array}[]{lll}\displaystyle\min_{x\in
X}&f_{1}(x)+\mathfrak{f}_{j}(\bar{x})+\langle\nabla\mathfrak{f}_{j}(\bar{x}),x-\bar{x}\rangle\\\
\mbox{s.t.}&c_{1}(x)+\mathfrak{c}_{\ell}(\bar{x})+\langle\nabla\mathfrak{c}_{\ell}(\bar{x}),x-\bar{x}\rangle\leq
0.\end{array}\right.$

for all $j\in A^{0}_{f}(\bar{x})$ and $\ell\in A^{0}_{c}(\bar{x})$. The result
follows from (Ackooij_Oliveira_2022, , Lem. 2.1). $\square$

###### Theorem 6.1

Consider model $\mathcal{M}$ in (15) and assume hypothesis (i) and (ii) in
Theorem 4.1 hold. Then every cluster point $\bar{x}$ of the sequence
$\\{x^{k}\\}_{k}$ computed by Algorithm 1 is M-critical. Furthermore, if
$\bar{x}$ is feasible and the pointwise Slater CQ holds, then $\bar{x}$ is
B-stationary.

Proof The result will follow from Theorem 4.1 once hypothesis (iii) is proven.
To this end, let $\\{x^{k_{i}}\\}_{i}$ be a subsequence of $\\{x^{k}\\}_{k}$
such that $\lim_{i}x^{k_{i}}=\bar{x}$. Observe that continuity of the convex
models $\mathcal{M}_{j\ell}$ and $\varepsilon>0$ fixed in the model’s
definition ensure the existence of a neighborhood $\mathcal{N}$ of $\bar{x}$
such that $A_{f}^{0}(\bar{x})\subset A_{f}^{\varepsilon}(y)$ and
$A_{c}^{0}(\bar{x})\subset A_{c}^{\varepsilon}(y)$ for all
$y\in\mathcal{N}\cap X$ (such set has nonempty relative interior as $\bar{x}$
belongs to the interior of $\mathcal{N}$ and $X$ is convex). Therefore, for
$i$ large enough,

$\displaystyle\mathcal{M}(y;x^{k_{i}})$ $\displaystyle=\displaystyle\min_{j\in
A_{f}^{\varepsilon}(x^{k_{i}})}\displaystyle\min_{\ell\in
A_{c}^{\varepsilon}(x^{k_{i}})}\mathcal{M}_{j\ell}(y;x^{k_{i}})\leq\displaystyle\min_{j\in
A_{f}^{0}(\bar{x})}\displaystyle\min_{\ell\in
A_{c}^{0}(\bar{x})}\mathcal{M}_{j\ell}(y;x^{k_{i}}).$

Let $x\in X$ be arbitrary, and $\\{y^{i}\\}_{i}\subseteq X$ such that
$y^{i}\to x$. By replacing $y$ with $y^{i}$ above and passing to the limit
superior with $i$ going to infinity we get, due to continuity of
$\mathcal{M}_{j\ell}$ in both arguments,

$\limsup_{i}\mathcal{M}(y^{i};x^{k_{i}})\leq\displaystyle\min_{j\in
A_{f}^{0}(\bar{x})}\displaystyle\min_{\ell\in
A_{c}^{0}(\bar{x})}\limsup_{i}\mathcal{M}_{j\ell}(y^{i};x^{k_{i}})=\displaystyle\min_{j\in
A_{f}^{0}(\bar{x})}\displaystyle\min_{\ell\in
A_{c}^{0}(\bar{x})}\mathcal{M}_{j\ell}(x;\bar{x})=\bar{\mathcal{M}}(x;\bar{x}).$

Noticing that $A_{H}(\bar{x})=A^{0}_{f}(\bar{x})\times A_{c}^{0}(\bar{x})$ in
Definition 2, and recalling item (ii) in Lemma 2, the result follows. Under
the extra assumption of $\bar{x}$ being feasible and holding Slater CQ, we
conclude from Lemma 4. $\square$

##  7 Optimization problems from classification applications

In this section, we consider problem (1) with functions given by

$f(x):=\sum_{j=1}^{N_{f}}G(\mathfrak{f}^{1}_{j}(x)+\mathfrak{f}_{j}^{2}(x)),\quad
c(x):=\sum_{j=1}^{N_{c}}R(\mathfrak{c}^{1}_{j}(x)+\mathfrak{c}^{2}_{j}(x)),$
(20)

with $\mathfrak{f}^{1}_{j},\mathfrak{c}^{1}_{j}:\mathbb{R}^{n}\to\mathbb{R}$
being convex functions,
$\mathfrak{f}^{2}_{j},\mathfrak{c}^{2}_{j}:\mathbb{R}^{n}\to\mathbb{R}$ being
concave, and $G,R:\mathbb{R}\to\mathbb{R}$ being convex and nondecreasing one-
dimensional functions. Problems of this form appear in machine learning
applications such as classification (RoysetWets_2022, , §6.I).

Given $\zeta_{j}\in\partial^{\tt C}\mathfrak{f}^{2}_{j}(x)$ and
$\xi_{j}\in\partial^{\tt C}\mathfrak{c}^{2}_{j}(x)$, we define the following
model

$\mathcal{M}_{a}(y;x):=\max\left\\{\begin{array}[]{lll}\sum_{j=1}^{N_{f}}G(\mathfrak{f}^{1}_{j}(y)+\mathfrak{f}_{j}^{2}(x)+\langle\zeta_{j},y-x\rangle)-\tau_{f}(x)\\\
\sum_{j=1}^{N_{c}}R(\mathfrak{c}^{1}_{j}(y)+\mathfrak{c}_{j}^{2}(x)+\langle\xi_{j},y-x\rangle)\end{array}\right\\},$
(21)

where $a\in A(x)\subseteq\left\\{(\zeta,\xi):\zeta_{{j}}\in\partial^{\tt
C}\mathfrak{f}^{2}_{j}(x)\;,\;\xi_{{j}}\in\partial^{\tt
C}\mathfrak{c}^{2}_{j}(x)\mbox{ for all }j\right\\}$. Thus, we consider the
model $\mathcal{M}(y;x)=\min_{a\in A(x)}\mathcal{M}_{a}(y;x)$ by assuming that
$1\leq|A(x)|\leq{\tt n}$ (cf. condition (8)). We see that Assumption 1(a)
holds. Note also that, as
$\mathfrak{f}^{1}_{j}(y)+\mathfrak{f}^{2}_{j}(y)\leq\mathfrak{f}^{1}_{j}(y)+\mathfrak{f}_{j}^{2}(x)+\langle\zeta_{j}(x),y-x\rangle$
and $G$ is nondecreasing (equivalently for $c$ and $R$), Assumption 1(b) holds
for $\bar{\mu}=0$.

###### Remark 3

As we clarified in Remark 2, given a cluster point $\bar{x}$, we can always
ensure (9) for the sequences of subgradients composing model (21). Considering
$\bar{\zeta}_{j}\in\partial^{\tt C}\mathfrak{f}^{2}_{j}(\bar{x})$ and
$\bar{\xi}_{j}\in\partial^{\tt C}\mathfrak{c}^{2}_{j}(\bar{x})$, the
corresponding cluster subgradient vectors for every $j$, we define
$A(\bar{x})$ to be the set of such cluster points provided by the subsequences
$\\{\zeta_{j}^{k_{i}}\\}_{i}$ and $\\{\xi_{j}^{k_{i}}\\}_{i}$:

$\bar{\mathcal{M}}(x;\bar{x})=\min_{a\in
A_{H}(\bar{x})}\max\left\\{\begin{array}[]{lll}\sum_{j=1}^{N_{f}}G(\mathfrak{f}^{1}_{j}(x)+\mathfrak{f}_{2}(\bar{x})+\langle\bar{\zeta}_{j},x-\bar{x}\rangle)-\tau_{f}(x)\\\
\sum_{j=1}^{N_{c}}R(\mathfrak{c}^{1}_{j}(x)+\mathfrak{c}_{2}(\bar{x})+\langle\bar{\xi}_{j},x-\bar{x}\rangle)\end{array}\right\\}.$

As in Remark 2, we can take
$A_{H}(\bar{x})=A(\bar{x})=\operatorname{LimOut}\left(A(x^{k_{i}})\right)$ for
the subsequence $\\{x^{k_{i}}\\}_{i}$ converging to $\bar{x}$.
$\blacktriangleleft$

###### Theorem 7.1

Consider Algorithm 1 furnished with model (21) applied to problem (1) with $f$
and $c$ given in (20). Assume that the sequence $\\{x^{k}\\}_{k}$ computed by
the algorithm is bounded. Then, any cluster point of $\\{x^{k}\\}_{k}$ is
M-critical.

Proof As in the proof of Theorem 5.1, we see that
$\lim_{i}\mathcal{M}(y^{i};x^{k_{i}})=\bar{\mathcal{M}}(y;\bar{x})$ for any
sequence $\\{y^{i}\\}_{i}\subseteq X$ such that $y^{i}\to y$. By relying on
Lemma 2 item (ii) we get hypothesis (iii) in Theorem 4.1. Hence, this last
theorem gives the stated result. $\square$

##  8 Numerical experiments

This section presents some preliminary numerical experiments to assess the
practical performance of our approach. We investigate two different (classes
of) stochastic problems fitting the setting in (5). A chance-constrained model
is considered for a gas network problem from LukMartHolgRene_2018 , and a
buffered probability optimization model for the design of a cantilever beam-
bar system studied in Song_2003  and Byun_Oliveira_Royset_2023 .

Our numerical experiments are performed in MATLAB R2022b on a computer Intel
Core i9 2.50GHz under Windows 11 Pro. The convex optimization solver called at
every iteration of our approach is the bundle algorithm described in Appendix
A. The algorithm’s parameters are set differently for each family of problems.

###  8.1 Chance-constrained optimization model

Let us consider that a random vector ${\boldsymbol{\omega}}$ belongs in a
certain probability space. Given a function $G$ depending on both
$x\in\mathbb{R}^{n}$ and $\boldsymbol{\omega}$, we say that $x$ is a success
for a given scenario $\omega$ if $G(x;\omega)\leq 0$ and it is a failure
otherwise. The uncertainty arising from this problem may be dealt with the
chance constraint

$\mathbb{E}[{\boldsymbol{1}}_{(0,+\infty)}(G(x;{\boldsymbol{\omega}}))]=\mathbb{P}[G(x;{\boldsymbol{\omega}})>0]\leq\alpha$

for a given $\alpha\in(0,1)$, which means that the probability of failure for
the decision variable $x$ is allowed to be, at most, the threshold $\alpha$.
Chance constraints arise in many real-life optimization problems, such as
electricity network expansion, energy management, mineral blending, structural
systems and other operating engineering systems
vanAckooij_Henrion_Moller_Zorgati_2011b ;
vanAckooij_Berge_Oliveira_Sagastizabal_2017 ; Prekopa_2003 ;
Rockafellar_Royset_2010 . The probabilistic constraint is traditionally faced
through the Sample Average Approximation (SAA) (see (SDR09, , §5)), which
leads to deterministic constraint as follows

$\sum_{j=1}^{N}p_{j}{\boldsymbol{1}}_{(0,+\infty)}[G(x;\omega^{j})]\leq\alpha,$
(22)

where $\\{\omega^{j}\\}_{j=1}^{N}$ is a set of $N$ realizations of the random
vector ${\boldsymbol{\omega}}$, and $p_{j}$ is the probability of scenario
$\omega^{j}$. It is clear that the above constraint does not fit (5) as
${\boldsymbol{1}}_{(0,+\infty)}$ is not even a continuous function.
Nevertheless, this issue will be overcome using a well-known nondecreasing
$\mathcal{C}^{2}$ approximation of ${\boldsymbol{1}}_{(0,+\infty)}$,
$\psi_{\theta}:\mathbb{R}\to(0,1)$ for $\theta>0$, given by

$\psi_{\theta}(x)=\frac{1}{1+\exp\left(-\frac{x}{\theta}\right)}.$

In the following example, we will take $\theta=0.1$ and a sample with
$N=10000$ scenarios.

#### Gas problem:

In this example, we consider an $n$-node connected tree (a graph with no
loops) representing a gas exit network as described in ClauHolRene_2016 ;
LukMartHolgRene_2018  (see Figure 1). The origin node (labeled as $0$) is the
one injecting the gas that flows through the pipes to the other nodes, where
the gas exits with a certain pressure. The exit load of each node $\ell$ is
random, and for each edge, a pressure drop is attributed. The goal is to find
a reliable upper bound for the pressure flowing out of each node, while
minimizing the sum of them. Let $\omega^{j}$ be a realization of the random
vector $\boldsymbol{\omega}$ representing the exit loads. A vector
$x\in\mathbb{R}^{n}$ is considered a success (for such a realization
$\omega^{j}$) if the inequalities

$G_{\ell}(x,{\omega^{j}}):=v_{\ell}({\omega^{j}})-x_{\ell}^{2}\leq 0$ (23)

hold for all $\ell=0,\dots,n-1$. Where
$v_{0}(\boldsymbol{\omega})=1+\max_{\ell=1,\dots,n-1}\\{h_{\ell}(\boldsymbol{\omega})\\}$,
$v_{\ell}(\boldsymbol{\omega})=v_{0}(\boldsymbol{\omega})-h_{\ell}(\boldsymbol{\omega})$
for $\ell=1,\dots,n-1$ and $h_{\ell}(\boldsymbol{\omega})$ is a function of
the random vector $\boldsymbol{\omega}$ (see ClauHolRene_2016 ;
LukMartHolgRene_2018 ).

Figure 1: The left-hand subfigure represents a graph with $n=4$ nodes for the
gas problem. The right-hand subfigure represents a graph with $n=12$ nodes.

The considered chance-constrained model is as follows

$\left\\{\begin{array}[]{lllll}\underset{x\geq
1}{\min}&\sum_{\ell=0}^{n-1}\;x_{\ell}\\\
\text{s.t.}&\displaystyle\sum_{j=1}^{N}\frac{1}{N}\max\left\\{\psi_{\theta}\big{(}G_{0}(x,\omega^{j})\big{)},\dots,\psi_{\theta}\big{(}G_{n-1}(x,\omega^{j})\big{)}\right\\}\leq\alpha,\end{array}\right.$

where the set of values $\\{h_{\ell}(\omega^{j})\\}_{j,\ell}$ have been
provided by LukMartHolgRene_2018 . Thus, the problem fits (5) as
$X=\\{x\in\mathbb{R}^{n}:x\geq 1\\}$ is an convex closed set. We solve it
employing the model described in (7) and compare our results against those
obtained from the DC approach of Javal_2021 . This reference approximates the
problem by a DC-constrained optimization model and applies the DC bundle
solver PBMDC2.

Table 1 reports the obtained results for two instances with $n=4$ and $n=12$
nodes. The following parameter values were set in Algorithm 1: $\kappa=0.3$,
$\lambda=0.1$, $\mu_{0}=2$ and $\texttt{Tol}=10^{-6}$.

$n$ | Iter | SS | CPU | $c(\bar{x})$ | $f(\bar{x})$
---|---|---|---|---|---
| Algorithm 1 | PBMDC2 |  | Algorithm 1 | PBMDC2 |  | Algorithm 1 | PBMDC2
$4$ | 3612 | 190 | 3600 | 113s | 53s | $-4.33\times 10^{-2}$ | $739.811$ | $738.797$
$12$ | 18185 | 956 | 18184 | 11min | 9min | $-1.26\times 10^{-4}$ | $3146.584$ | $3142.05$
Table 1: Computed results for the Gas Problem for $n=4$ and $n=12$ compared to
those provided by Algorithm PBMDC2 from Javal_2021 . The Iter and SS values
show the number of iterations and serious steps computed. Here, $\bar{x}$ is
the last iterate computed by the algorithms.

We have run PBMDC2 with the same parameters described in Javal_2021 . We
refrain from comparing our algorithm against the approach of
LukMartHolgRene_2018  since it has been shown in Javal_2021  that PBMDC2
performs better for this problem. Although our method is more general, Table 1
shows that it provides good results compared to a more specific algorithm.

###  8.2 Buffered probability optimization model

We consider a buffered failure probability constraint on an integrable random
variable ${\boldsymbol{\xi}}$ and recall that, for $\alpha\in(0,1)$ one has
Rockafellar_Royset_2010 :

buffered probability of ${\boldsymbol{\xi}}>0$ is no greater than $1-\alpha$
if and only if ${\tt AVaR}_{\alpha}({\boldsymbol{\xi}})\leq 0$,

where ${\tt AVaR}_{\alpha}({\boldsymbol{\xi}})$ is the _Average value-at-Risk_
of ${\boldsymbol{\xi}}$, given by

$\min_{t\in\mathbb{R}}\;t+\frac{1}{1-\alpha}\mathbb{E}[\max\left\\{0,{\boldsymbol{\xi}}-t\right\\}].$

For further details about buffered failure probabilities, we refer to
(RoysetWets_2022, , § 3.E) and Mafusalov_Uryasev_2018 .

Let us consider that
${\boldsymbol{\xi}}=\max_{\ell=1,\dots,m}\mathfrak{c}_{\ell}(y;{\boldsymbol{\omega}})$,
with $\mathfrak{c}_{j}(\cdot;{\boldsymbol{\omega}})$ being a weakly-concave
function. After sampling according to the distribution of
${\boldsymbol{\omega}}$, as we did in (22), we obtain the following expression
for the constraint ${\tt AVaR}_{\alpha}({\boldsymbol{\xi}})\leq 0$

$c(y,t)=-t\frac{\alpha}{1-\alpha}+\sum_{j=1}^{N}\frac{p_{j}}{1-\alpha}\;\max\\{t,\mathfrak{c}_{1}(y,\omega^{j}),\dots,\mathfrak{c}_{m}(y,\omega^{j})\\}\leq
0.$

Denoting $c_{1}(y,t)=-t\frac{\alpha}{1-\alpha}$ and
$c_{2}(y,t)=\sum_{j=1}^{N}\frac{p_{j}}{1-\alpha}\;\max\\{t,\mathfrak{c}_{1}(y,\omega^{j}),\dots,\mathfrak{c}_{m}(y,\omega^{j})\\}$,
we end up in a suitable expression for problem (5), with no need of a second
layer of approximation. It is worth noting that ${\tt
AVaR}_{\alpha}({\boldsymbol{\xi}})$ operator possesses a relevant
characteristic of being level-bounded on $X\times\mathbb{R}$. This ensures
that the sequence generated by the algorithm is bounded.

#### Design of cantilever beam-bar system:

This example investigates an optimal design of a cantilever beam-bar system
illustrated in Figure 2. The system consists of an ideally plastic cantilever
beam of moment capacity $\boldsymbol{\omega}_{M}+y_{M}$ and length $2L$, with
$L=5$, which is propped by an ideally rigid-brittle bar of strength
$\boldsymbol{\omega}_{T}+y_{T}$. The structure is subjected to load
$\boldsymbol{\omega}_{P}$ that is applied in the middle of the bar. The moment
capacity, strength, and load are assumed to be Gaussian random variables:
$\boldsymbol{\omega}_{M}+y_{M}\sim N(y_{M},300)$,
$\boldsymbol{\omega}_{T}+y_{T}\sim N(y_{T},20)$, and
$\boldsymbol{\omega}_{P}\sim N(150,30)$. In this example, we optimize the two
design variables $y_{M}\in[500,1500]$ and $y_{T}\in[50,150]$, representing the
expected values of the moment capacity and the strength, respectively. The
cost function is $2y_{M}+y_{T}$. Then, the system can be described using the
following limit-state functions determining whether a component fails or
survives:

$\displaystyle g_{1}(y,\boldsymbol{\omega}):=$
$\displaystyle-(y_{T}+\boldsymbol{\omega}_{T})+\frac{5}{16}\boldsymbol{\omega}_{P},$
$\displaystyle g_{2}(y,\boldsymbol{\omega}):=$
$\displaystyle\vphantom{\frac{1}{1}}-(y_{M}+\boldsymbol{\omega}_{M})+L\boldsymbol{\omega}_{P},$
$\displaystyle g_{3}(y,\boldsymbol{\omega}):=$
$\displaystyle-(y_{M}+\boldsymbol{\omega}_{M})+\frac{3L}{8}\boldsymbol{\omega}_{P},$
$\displaystyle g_{4}(y,\boldsymbol{\omega}):=$
$\displaystyle-(y_{M}+\boldsymbol{\omega}_{M})+\frac{L}{3}\boldsymbol{\omega}_{P},$
$\displaystyle g_{5}(y,\boldsymbol{\omega}):=$
$\displaystyle\vphantom{\frac{1}{1}}-(y_{M}+\boldsymbol{\omega}_{M})-2L(y_{T}+\boldsymbol{\omega}_{T})+L\boldsymbol{\omega}_{P}.$

Figure 2: Example cantilever beam-bar system. Figure extracted from
Byun_Oliveira_Royset_2023 .

The system’s failure can be thus represented as
$\max\\{G_{1}(y,{\omega^{j}}),G_{2}(y,{\omega^{j}}),G_{3}(y,{\omega^{j}})\\}>0$,
with

$\displaystyle G_{1}(y,\boldsymbol{\omega}):=$
$\displaystyle\min\\{g_{1}(y,\boldsymbol{\omega}),g_{2}(y,\boldsymbol{\omega})\\},$
$\displaystyle G_{2}(y,\boldsymbol{\omega}):=$
$\displaystyle\min\\{g_{3}(y,\boldsymbol{\omega}),g_{4}(y,\boldsymbol{\omega})\\},$
$\displaystyle G_{3}(y,\boldsymbol{\omega}):=$
$\displaystyle\min\\{g_{3}(y,\boldsymbol{\omega}),g_{5}(y,\boldsymbol{\omega})\\}.$

Our buffered probability model for the problem reads as

$\left\\{\begin{array}[]{llll}\underset{(y,t)\in X}{\min}&2y_{M}+y_{T}\\\
\text{s.t.}&\displaystyle-t\frac{\alpha}{1-\alpha}+\sum_{j=1}^{N}\frac{1}{N(1-\alpha)}\max\left\\{G_{1}(y,\omega^{j}),G_{2}(y,\omega^{j}),G_{3}(y,\omega^{j}),t\right\\}\leq
0,\end{array}\right.$

where $\omega^{j}$ is a realization of the random vector
$\boldsymbol{\omega}=(\boldsymbol{\omega}_{M},\boldsymbol{\omega}_{T},\boldsymbol{\omega}_{P})$,
$y=(y_{M},y_{T})$ and

$X=\left\\{(y,t)\in\mathbb{R}^{3}:\begin{array}[]{l}500\leq y_{M}\leq 1500\\\
50\leq y_{T}\leq 150\end{array}\right\\}.$

We consider the model described in (7) to solve problems such as (5). For this
problem, we have considered a sample with $N=100000$ scenarios and
$\alpha=0.999$. The algorithm’s parameters were set as $\kappa=0.3$,
$\lambda=0.1$, $\mu_{0}=\kappa$ and $\texttt{Tol}=10^{-6}$. The algorithm
stops after 35 seconds and 183 iterations (all of them Serious Steps). The
computed feasible point $(1288,150)$ has a cost of $2727$, which is slightly
smaller (under $0.1$% of relative error) than the reference value $2729$ for
the feasible point $(1289,150)$, issued by a grid search on the set $X$ with
$1000\times 100$ points. Such grid search took about 240 seconds of
computation.

###### Acknowledgements.

The second author acknowledges financial support from the Gaspard-Monge
program for Optimization and Operations Research (PGMO) project “SOLEM -
Scalable Optimization for Learning and Energy Management”. The third author
acknowledges support from the Office of Naval Research.

## References

  * (1) J. V. Burke and M. C. Ferris.  A Gauss—Newton method for convex composite optimization.  Mathematical Programming, 71(2):179–194, 1995.
  * (2) C. Sagastizábal.  Composite proximal bundle method.  Mathematical Programming, 140(1):189–233, 2013.
  * (3) A. S. Lewis and S. J. Wright.  A proximal method for composite minimization.  Mathematical Programming, 158(1):501–546, 2016.
  * (4) Y. Cui, J-S Pang, and B. Sen.  Composite difference-max programs for modern statistical estimation problems.  SIAM Journal on Optimization, 28(4):3344–3374, 2018.
  * (5) J. Liu, Y. Cui, J-S Pang, and S. Sen.  Two-stage stochastic programming with linearly bi-parameterized quadratic recourse.  SIAM Journal on Optimization, 30(3):2530–2558, 2020.
  * (6) Y. Cui, J. Liu, and J-S Pang.  Nonconvex and nonsmooth approaches for affine chance-constrained stochastic programs.  Set-Valued and Variational Analysis, 30(3):1149–1211, 2022.
  * (7) F.J. Aragón-Artacho, P. Pérez-Aros, and D. Torregrosa-Belén.  The boosted double-proximal subgradient algorithm for nonconvex optimization, 2023.
  * (8) J. V. Burke.  Descent methods for composite nondifferentiable optimization problems.  Mathematical Programming, 33(3):260–279, 1985.
  * (9) C. Sagastizábal.  Divide to conquer: Decomposition methods for energy optimization.  Mathematical Programming, 134(1):187–222, 2012.
  * (10) L. T. H. An and P. D. Tao.  The DC (difference of convex functions) programming and DCA revisited with DC models of real world nonconvex optimization problems.  Annals of Operations Research, 133(1):23–46, 2005.
  * (11) W. de Oliveira.  Proximal bundle methods for nonsmooth dc programming.  Journal of Global Optimization, 75(2):523–563, 2019.
  * (12) W. van Ackooij and W. de Oliveira.  Some brief observations in minimizing the sum of locally lipschitzian functions.  Optimization Letters, 14(3):509–520, 2020.
  * (13) J. O. Royset.  Consistent approximations in composite optimization.  Mathematical Programming, 201(1):339–372, 2023.
  * (14) Y. Cui and J-S Pang.  Modern Nonconvex Nondifferentiable Optimization.  SIAM, 2022.
  * (15) P. Apkarian, D. Noll, and A. Rondepierre.  Mixed $H_{2}/H_{\infty}$ Control via Nonsmooth Optimization.  SIAM Journal on Control and Optimization, 47(3):1516–1546, 2008\.
  * (16) H. A. Le Thi, T. Pham Dinh, and H. V. Ngai.  Exact penalty and error bounds in dc programming.  Journal of Global Optimization, 52(3):509–535, 2012.
  * (17) A. Strekalovsky and I. Minarchenko.  On local search in d.c. optimization.  In 2017 Constructive Nonsmooth Analysis and Related Topics (dedicated to the memory of V.F. Demyanov) (CNSA), pages 1–4, May 2017.
  * (18) J-S Pang, M. Razaviyayn, and A. Alvarado.  Computing b-stationary points of nonsmooth dc programs.  Mathematics of Operations Research, 42(1):95–118, 2017.
  * (19) O. Montonen and K. Joki.  Bundle-based descent method for nonsmooth multiobjective DC optimization with inequality constraints.  Journal of Global Optimization, 72(3):403–429, 2018.
  * (20) Y. Cui, J-S Pang, and B. Sen.  Composite difference-max programs for modern statistical estimation problems.  SIAM Journal on Optimization, 28(4):3344–3374, 2018.
  * (21) W. van Ackooij, S. Demassey, P. Javal, H. Morais, W. de Oliveira, and B. Swaminathan.  A bundle method for nonsmooth DC programming with application to chance-constrained problems.  Computational Optimization and Applications, 78(2):451–490, 2021\.
  * (22) K. Syrtseva, W. de Oliveira, S. Demassey, and W. van Ackooij.  Minimizing the difference of convex and weakly convex functions via bundle method.  Pacific Journal of Optimization, 2024.
  * (23) F.J. Aragón-Artacho, B.S. Mordukhovich, and P. Pérez-Aros.  Coderivative-based semi-newton method in nonsmooth difference programming, 2023.
  * (24) J-E Byun, W. de Oliveira, and J. O. Royset.  S-BORM: Reliability-based optimization of general systems using buffered optimization and reliability method.  Reliability Engineering & System Safety, 236:109314, 2023.
  * (25) W. van Ackooij and W. de Oliveira.  Non-smooth dc-constrained optimization: constraint qualification and minimizing methodologies.  Optimization Methods and Software, 34(4):890–920, 2019.
  * (26) C. Sagastizábal and M. Solodov.  An infeasible bundle method for nonsmooth convex constrained optimization without a penalty function or a filter.  SIAM Journal on Optimization, 16(1):146–169, 2005.
  * (27) F.H. Clarke.  Optimisation and Nonsmooth Analysis.  Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 1987.
  * (28) J.O. Royset and R.J-B Wets.  An optimization primer.  Springer Series in Operations Research and Financial Engineering. Springer Cham, 1 edition, 2022.
  * (29) J.L. Higle and S. Sen.  On the convergence of algorithms with implications for stochastic and nondifferentiable optimization.  Mathematics of operations research, 17(1):112–131, 1992.
  * (30) R. T. Rockafellar and R. J-B Wets.  Variational Analysis, volume 317 of Grundlehren der mathematischen Wissenschaften.  Springer Verlag Berlin, 3rd edition, 2009.
  * (31) K. Syrtseva, W. de Oliveira, S. Demassey, H. Morais, P. Javal, and B. Swaminathan.  Difference-of-convex approach to chance-constrained optimal power flow modelling the dso power modulation lever for distribution networks.  Sustainable Energy, Grids and Networks, 36, 2023.
  * (32) R. T. Rockafellar and R. J-B Wets.  A Lagrangian finite generation technique for solving linear-quadratic problems in stochastic programming.  Mathematical Programming Study, 28:63–93, 1986.
  * (33) J. S. Pang.  Partially B-regular optimization and equilibrium problems.  Mathematics of Operations Research, 32(3):687–699, 2007.
  * (34) W. van Ackooij and W. de Oliveira.  Addendum to the paper ‘nonsmooth DC-constrained optimization: constraint qualification and minimizing methodologies’.  Optimization Methods and Software, pages 2241–2250, 2022.
  * (35) R. Correa and C. Lemaréchal.  Convergence of some algorithms for convex minimization.  Mathematical Programming, 62(2):261–275, 1993.
  * (36) L. Adam, M. Branda, H. Heitsch, and R. Henrion.  Solving joint chance constrained problems using regularization and benders’ decomposition.  Annals of Operations Research, 292:683–709, 2018.
  * (37) J. Song and A. Der Kiureghian.  Bounds on system reliability by linear programming.  Journal of Engineering Mechanics, 129(6):627–636, 2003.
  * (38) W. van Ackooij, R. Henrion, A. Möller, and R. Zorgati.  Joint chance constrained programming for hydro reservoir management.  Optimization and Engineering, 15:509–531, 2014.
  * (39) W. van Ackooij, V. Berge, W. de Oliveira, and C. Sagastizábal.  Probabilistic optimization via approximate p-efficient points and bundle methods.  Computers & Operations Research, 77:177 – 193, 2017.
  * (40) A. Prékopa.  Probabilistic programming.  In A. Ruszczyński and A. Shapiro, editors, Stochastic Programming, volume 10 of Handbooks in Operations Research and Management Science, pages 267–351. Elsevier, Amsterdam, 2003.
  * (41) R.T. Rockafellar and J.O. Royset.  On buffered failure probability in design and optimization of structures.  Reliability Engineering & System Safety, 95(5):499–510, 2010.
  * (42) A. Shapiro, D. Dentcheva, and A. Ruszczyński.  Lectures on Stochastic Programming: Modeling and Theory.  MPS-SIAM Series on Optimization. SIAM - Society for Industrial and Applied Mathematics and Mathematical Programming Society, Philadelphia, 2009.
  * (43) C. Gotzes, H. Heitsch, and R. Henrion.  On the quantification of nomination feasibility in stationary gas networks with random load.  Mathematical Methods of Operational Research, 84:427–457, 2016\.
  * (44) A. Mafusalov and S. Uryasev.  Buffered probability of exceedance: Mathematical properties and optimization.  SIAM Journal on Optimization, 28(2):1077–1103, 2018.
  * (45) Welington de Oliveira and Mikhail Solodov.  Bundle Methods for Inexact Data, pages 417–459.  Springer International Publishing, Cham, 2020.
  * (46) A. Frangioni.  Standard Bundle Methods: Untrusted Models and Duality, pages 61–116.  Springer International Publishing, Cham, 2020.

##  Appendix A On the (inexact) solution of subproblems: a prox-form of bundle
method

In this section, we rely on Correa_Lemarechal_1993  and show how a bundle-like
algorithm can be deployed for computing an $\epsilon^{k}$-solution $y^{k}$ of
(3) for the class of models $\mathcal{M}$ holding Assumption 1. Note that, in
general, the chosen model need not be convex. Following the idea in Remark 1,
we propose an algorithm to solve every model $\mathcal{M}_{a}(\cdot;x^{k})$,
which happens to be convex.

In what follows, let $k$ be a fixed iteration of Algorithm 1, $x^{k}$ be the
current point, and $\mathcal{M}_{a}(\cdot;x^{k})$ for a fixed $a\in A(x^{k})$.
In the following approach, iterations will be denoted by $m$: our goal is to
generate a sequence of _inner iterates_ $\\{y^{k,m}_{a}\\}_{m}\subset X$ with
the following properties:

  1. (a)

if $x^{k}$ solves (4), then $\lim_{m\to\infty}y^{k,m}_{a}=x^{k}$;

  2. (b)

if $x^{k}$ is not a solution to (4), then $\\{y^{k,m}_{a}\\}_{m}$ is finite
and its last term, denoted by $y^{k}_{a}$, is an $\epsilon^{k}_{a}$-solution
of (4), with $\epsilon^{k}_{a}$ satisfying
$\epsilon^{k}_{a}\leq\frac{\lambda}{2}\textstyle{\|}y^{k}_{a}-x^{k}\textstyle{\|}^{2}.$

To this end, at every inner iteration $m$, we approximate a convex function
$\mathcal{M}_{a}(\cdot;x^{k})$ with a cutting-plane model
$\check{\mathcal{M}}^{k,m}_{a}\leq\mathcal{M}_{a}(\cdot;x^{k})$ and define
$y^{k,m+1}_{a}$ as

$y^{k,m+1}_{a}:=\operatorname{argmin}_{x\in
X}\;\check{\mathcal{M}}^{k,m}_{a}(x)+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}.$
(24)

As we will precise below, $\check{\mathcal{M}}^{k,m}_{a}$ is given by at most
$m$ pieces. Thus, if $X$ is polyhedral, then the above problem is a strictly
convex quadratic programming problem, for which efficient (open source and
commercial) solvers exist. As standard in bundle methods (see for instance
Oliveira_Solodov_2020 ; Frangioni_2020 ) the number of cuts in the model
$\check{\mathcal{M}}^{k,m}_{a}$ can be bounded. Indeed, only two pieces are
enough to have a convergent algorithm. In more general terms, any rule to
update $\check{\mathcal{M}}^{k,m}_{a}$ to $\check{\mathcal{M}}^{k,m+1}_{a}$ is
a valid one provided that the following inequalities are satisfied:

$\begin{matrix}\hfill\check{\mathcal{M}}^{k,m+1}_{a}\leq&\mathcal{M}_{a}(\cdot;x^{k}),\\\
\check{\mathcal{M}}^{k,m}_{a}(y^{k,m+1}_{a})+\mu^{k}\langle
x^{k}-y^{k,m+1}_{a},\cdot-y^{k,m+1}_{a}\rangle\leq&\check{\mathcal{M}}^{k,m+1}_{a},\\\
\hfill\mathcal{M}_{a}(y^{k,m+1}_{a},x^{k})+\langle s^{k,m}_{a},\cdot-
y_{a}^{k,m+1}\rangle\leq&\check{\mathcal{M}}^{k,m+1}_{a},\end{matrix}$ (25)

where $s^{k,m}_{a}\in\partial\mathcal{M}_{a}(\cdot;x^{k})(y^{k,m}_{a})$ is an
arbitrary subgradient.

######  Theorem A.1 ((Correa_Lemarechal_1993, , Proposition 4.3))

Consider a sequence of convex models $\\{\check{\mathcal{M}}^{k,m}_{a}\\}_{m}$
such that conditions (25) hold and sequence $\\{y^{k,m}_{a}\\}_{m}$ given by
(24). Then,

$\lim_{m\to\infty}\mathcal{M}_{a}(y^{k,m+1}_{a};x^{k})-\check{\mathcal{M}}_{a}^{k,m}(y^{k,m+1}_{a})=0\quad\mbox{and}\quad\lim_{m\to\infty}y^{k,m}_{a}=\bar{y}^{k}_{a},$

where $\bar{y}^{k}_{a}$ is the unique solution to (4).

Given these ingredients, we now present the following algorithm for inexactly
solving (3).

Algorithm 2 prox-form of bundle method

1:Let $x^{k},\,\mu^{k},\,\kappa,\,\lambda$, and Tol be provided by Algorithm 1
at iteration $k$

2:Set $y^{k,0}_{a}=x^{k}$ and define
$\check{\mathcal{M}}^{k,0}_{a}:=\mathcal{M}_{a}(y_{a}^{k,0};x^{k})+\langle
s_{a}^{k,0},\cdot-y_{a}^{k,0}\rangle$, with
$s_{a}^{k,0}\in\partial\mathcal{M}_{a}(\cdot;x^{k})(y_{a}^{k,0})$

3:

4:for $m=0,1,2,\dots$ do

5: Let $y^{k,m+1}_{a}$ be the unique solution of  $\min_{x\in
X}\;\check{\mathcal{M}}_{a}^{k,m}(x)+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}$

6: if
$\mathcal{M}_{a}(x^{k};x^{k})-\check{\mathcal{M}}_{a}^{k,m}(y_{a}^{k,m+1})\leq\texttt{Tol}$
then

7: Stop and return $y^{k}_{a}:=x^{k}$

8: end if

9: if
$\mathcal{M}_{a}(y_{a}^{k,m+1};x^{k})-\check{\mathcal{M}}_{a}^{k,m}(y_{a}^{k,m+1})\leq\frac{\lambda}{2}\textstyle{\|}y_{a}^{k,m+1}-x^{k}\textstyle{\|}^{2}$
then

10: Stop and return $y_{a}^{k}:=y_{a}^{k,m+1}$

11: end if

12:

13: Update the model $\check{\mathcal{M}}_{a}^{k,m}$ to
$\check{\mathcal{M}}_{a}^{k,m+1}$ in order to satisfy (25)

14:end for

Let $s_{a}^{k,j}\in\partial\mathcal{M}_{a}(\cdot;x^{k})(y_{a}^{k,j})$,
$j=1,\ldots,m$. As standard in bundle methods, a practical rule to update the
model is as follows:

$\check{\mathcal{M}}_{a}^{k,m+1}(\cdot):=\max_{j\in
J_{m+1}}\\{\mathcal{M}_{a}(y^{k,j};x^{k})+\langle s_{a}^{k,j},\cdot-
y_{a}^{k,j}\rangle\\},$

where $J_{m+1}:=\\{{J^{act}_{m},m+1}\\}$, and $J^{act}_{m}\subset
J_{m}\subset\\{0,1,\ldots,m\\}$ is the index set of active constraints of the
following reformulation of problem (3):

$\left\\{\begin{array}[]{lll}\displaystyle\min_{y\in
X,r\in\mathbb{R}}&r+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}\\\
\mbox{s.t.}&\mathcal{M}_{a}(y_{a}^{k,j};x^{k})+\langle
s_{a}^{k,j},x-y_{a}^{k,j}\rangle\leq r&\forall\,j\in J_{m}.\end{array}\right.$

It is not difficult to show that this rule satisfies (25).

###### Theorem A.2

Consider Algorithm 2 with $\texttt{Tol}=0$ applied to problem (4).

  1. (i)

If the algorithm loops indefinitely, or it stops at line 7, then $x^{k}$ is a
solution to (4).

  2. (ii)

If the algorithm halts at line 10, then $y^{k}_{a}:=y_{a}^{k,m+1}$ is an
$\epsilon^{k}_{a}$-solution of (4), with
$\epsilon_{a}^{k}:=\mathcal{M}_{a}(y_{a}^{k,m+1};x^{k})-\check{\mathcal{M}}_{a}^{k,m}(y_{a}^{k,m+1})\geq
0$ satisfying
$\epsilon^{k}_{a}\leq\frac{\lambda}{2}\textstyle{\|}y^{k}_{a}-x^{k}\textstyle{\|}^{2}$.

  3. (iii)

Suppose that $x^{k}$ does not solve (4). Then Algorithm 2 stops after finitely
many iterations.

Proof. To prove item (i), note that if the algorithm stops at line 7, then

$\displaystyle\check{\mathcal{M}}_{a}^{k,m}(y_{a}^{k,m+1})+{\frac{\mu^{k}}{2}}\textstyle{\|}y_{a}^{k,m+1}-x^{k}\textstyle{\|}^{2}$
$\displaystyle=\min_{y\in
X}\left\\{\check{\mathcal{M}}_{a}^{k,m}(y)+\frac{\mu^{k}}{2}\textstyle{\|}y-x^{k}\textstyle{\|}^{2}\right\\}$
$\displaystyle\leq\min_{y\in
X}\left\\{\mathcal{M}_{a}(y;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y-x^{k}\textstyle{\|}^{2}\right\\}$
$\displaystyle\leq\mathcal{M}_{a}(x^{k};x^{k})\leq\check{\mathcal{M}}_{a}^{k,m}(y_{a}^{k,m+1})+\texttt{Tol}.$

In this case, $\texttt{Tol}=0$ implies that $y_{a}^{k,m+1}=x^{k}$ solves (4).
On the other hand, if the algorithm loops indefinitely, then

$\frac{\lambda}{2}\textstyle{\|}y^{k,m+1}_{a}-x^{k}\textstyle{\|}^{2}<\mathcal{M}_{a}(y_{a}^{k,m+1};x^{k})-\check{\mathcal{M}}_{a}^{k,m}(y_{a}^{k,m+1})\quad\forall\,m.$

Theorem A.1 then ensures that $\lim_{m\to\infty}y^{k,m}_{a}=x^{k}$ solves (4).

To show item (ii), note that by setting
$\epsilon_{a}^{k}:=\mathcal{M}_{a}(y_{a}^{k,m+1};x^{k})-\check{\mathcal{M}}_{a}^{k,m}(y^{k,m+1})\geq
0$ we get

$\displaystyle\mathcal{M}_{a}(y^{k,m+1}_{a};x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y_{a}^{k,m+1}-x^{k}\textstyle{\|}^{2}$
$\displaystyle=\check{\mathcal{M}}^{k,m}_{a}(y_{a}^{k,m+1})+{\frac{\mu^{k}}{2}}\textstyle{\|}y_{a}^{k,m+1}-x^{k}\textstyle{\|}^{2}+\epsilon^{k}_{a}$
$\displaystyle=\min_{x\in
X}\left\\{\check{\mathcal{M}}_{a}^{k,m}(x)+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}\right\\}+\epsilon^{k}_{a}$
$\displaystyle\leq\min_{x\in
X}\left\\{\mathcal{M}_{a}(x;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}\right\\}+\epsilon_{a}^{k},$

i.e., $y^{k,m+1}_{a}$ is an $\epsilon_{a}^{k}$-solution to (4). By assumption
the algorithm stops at line 10, then
$\epsilon^{k}_{a}\leq\frac{\lambda}{2}\textstyle{\|}y^{k,m+1}_{a}-x^{k}\textstyle{\|}^{2}$
and the result follows.

Finally, to prove item (iii), see that if $x^{k}$ does not solve (4), then
Theorem 2 ensures that for every $\epsilon>0$ small enough, there exists
$m_{0}\in\mathbb{N}$ such that

$\mathcal{M}_{a}(y^{k,m+1}_{a};x^{k})-\check{\mathcal{M}}^{k,m}_{a}(y_{a}^{k,m+1})<\epsilon\quad\mbox{and}\quad\epsilon<\frac{\lambda}{2}\textstyle{\|}y^{k,m+1}_{a}-x^{k}\textstyle{\|}^{2}\quad\mbox{for
all }\;m\geq m_{0},$

where $\lambda$ is the one in Algorithm 2. Hence, the stopping test at line 9
is triggered and the algorithm stops after finitely many steps. $\square$

###### Remark 4

We need to apply Algorithm 2 for all $a\in A(x^{k})$ in order to compute an
$\epsilon^{k}_{a}$-solution to every problem (4). Once all these iterations
are done, we compute $y^{k}$ as

$y^{k}\in\operatorname{argmin}_{y\in
Y^{k}}\mathcal{M}(y;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y-x^{k}\textstyle{\|}^{2},$

where $Y^{k}=\\{y^{k}_{a}:a\in A(x^{k})\\}$ is a finite set. We can easily
prove that $y^{k}$ is such that

$\mathcal{M}(y^{k};x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}\leq\mathcal{M}(y^{k}_{a};x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y^{k}_{a}-x^{k}\textstyle{\|}^{2}\leq\mathcal{M}_{a}(y^{k}_{a};x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y^{k}_{a}-x^{k}\textstyle{\|}^{2}\leq\mathcal{M}_{a}(x;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}+\epsilon^{k}_{a}$

for all $a\in A(x^{k})$ and $x\in X$. Then, minimizing over $x\in X$ and $a\in
A(x^{k})$ the right-hand side of the above chain of inequalities, we deduce

$\mathcal{M}(y^{k};x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}\leq\min_{x\in
X}\mathcal{M}(x;x^{k})+\frac{\mu^{k}}{2}\textstyle{\|}x-x^{k}\textstyle{\|}^{2}+\epsilon^{k},$

where $\epsilon^{k}=\min_{a\in A(x^{k})}\epsilon^{k}_{a}$. As
$y^{k}=y^{k}_{a^{*}}$ for a given $a^{*}\in A(x^{k})$, we conclude that
$y^{k}$ is an $\epsilon^{k}$-solution to problem (3) with

$\epsilon^{k}\leq\epsilon^{k}_{a^{*}}\leq\frac{\lambda}{2}\textstyle{\|}y^{k}-x^{k}\textstyle{\|}^{2}.$

$\blacktriangleleft$