# Mixed-integer linearity in nonlinear optimization:
a trust region approach

Alberto De Marchi

###### Abstract

Bringing together nonlinear optimization with mixed-integer linear constraints
enables versatile modeling, but poses significant computational challenges. We
investigate a method to solve these problems based on sequential mixed-integer
linearization with trust region safeguard, computing feasible iterates via
calls to a generic mixed-integer linear solver. Convergence to critical,
possibly suboptimal, feasible points is established for arbitrary starting
points. Finally, we present numerical applications in nonsmooth optimal
control and optimal network design and operation.University of the Bundeswehr
Munich, Department of Aerospace Engineering, Institute of Applied Mathematics
and Scientific Computing, 85577 Neubiberg, Germany.email
<EMAIL_ADDRESS>orcid 0000-0002-3545-6898.

Keywords. Mixed-integer programming, nonlinear programming, successive
linearization schemes, trust region methods.

AMS MSC. 65K05, 90C06, 90C11, 90C30.

##  1 Introduction

The operation of complex, integrated process systems demands efficient use of
resources whilst imposing tight safety constraints. Mixed-integer optimization
provides a powerful and flexible mathematical template for modeling many tasks
that involve discrete and continuous variables. An equivalent framework is
that of generalized disjunctive programming (GDP), which allows to easily
capture both algebraic relations and logical propositions; see [9]. For
transforming GDP into mixed-integer nonlinear programming (MINLP), Boolean
variables can be replaced by binary ones and logical specifications can be
converted, possibly in relaxed form, into linear algebraic constraints.

In this work we consider optimization problems with a nonlinear objective and
mixed-integer linear constraints

$\displaystyle\operatorname*{minimize}\quad f(x)\qquad\operatorname{over}\quad
x\in\mathcal{X}:=\left\\{x\in\mathbb{R}^{n}\,\middle|\,A_{\textup{e}}x=b_{\textup{e}},~{}x_{\textup{l}}\leq
x\leq x_{\textup{u}},~{}x_{i}\in\mathbb{Z}\quad\forall
i\in\mathcal{I}\right\\}$ (P)

with smooth objective function $f\colon\mathcal{X}\to\mathbb{R}$ and nonempty
closed feasible set $\mathcal{X}\subseteq\mathbb{R}^{n}$ described by
polyhedral and integrality constraints defined by some matrix
$A_{\textup{e}}$, vectors $b_{\textup{e}}$, $x_{\textup{l}}$,
$x_{\textup{u}}$, and index set $\mathcal{I}\subseteq\\{1,\ldots,n\\}$. We
intend to exploit the tractacle, although difficult, structure of
$\mathcal{X}$ in our numerical method, in particular avoiding quadratic terms
that may arise from local models or regularizations and instead relying on
mixed-integer linear programs (MILP) as subproblems. To do so, we make use of
_sharp_ norms, namely $\ell_{p}$ with $p\in\\{1,\infty\\}$, instead of the
classical Euclidean $\ell_{2}$. Then, a practical assumption is that MILPs
involving $\mathcal{X}$ can be efficiently solved, leveraging mature LP and
MILP technology. Moreover, we should mention that more general problems having
also nonlinear constraints could be tackled as well, by integrating the
approach for solving (P) within sequential partially-constrained optimization
schemes, such as augmented Lagrangian and interior point methods.

Partitioning the decision variables $x$ into real-valued $u$ and integer-
valued $z$, we consider hereafter the following blanket assumptions on (P):

  1. (i)

the feasible set $\mathcal{X}$ has mixed-integer linear structure;

  2. (ii)

the objective function $f$ is linear with respect to $z$;

  3. (iii)

all feasible values for $z$ lie in a bounded set.

Assumption (ii) of integer-linearity for $f$ poses some modelling restrictions
but avoids the need for derivative approximations or exotic notions thereof;
in line with [13], and in contrast to [7], we assume the objective $f$ to be
of the form $(u,z)\mapsto f_{1}(u)+\langle f_{2},z\rangle$, with $f_{2}$ being
a vector of suitable size. Notice that for the more general case with
additional nonlinear constraints Assumption (ii) would not be restrictive. By
requiring the integer-valued set to be bounded, Assumption (iii) is not only
reasonable and often valid in practice, but it also greatly simplifies the
presentation, possibly with only minor algorithmic consequences. Furthermore,
we investigate (P) with the implicit assumption that there are both real- and
integer-valued variables as well as some couplings thereof in $\mathcal{X}$
and nonlinearities in $f$, for otherwise off-the-shelf methods could be
readily adopted.

Focusing on (P) without any convexity assumptions, we will not demand _global_
optimality but seek instead an affordable local “solution”. The underlying
optimality notion is defined by localization in a trust region. We analyze the
similarities of projected gradient-type methods with trust region
linearization schemes, highlighting how sharp norms may lead to undesirable
properties for the former, making the latter an interesting alternative. After
detailing an iterative numerical scheme, we investigate its global convergence
properties, in the sense of arbitrary initial points. The algorithm generates
feasible iterates with (nonmonotonically) decreasing objective, and relies on
averaged merit values to monitor and enforce convergence.

#### Related Work

Based on the work of [4], sequential linearization methods appear to heavily
exploit convexity of $\mathcal{X}$, but the approach has been applied to
nonconvex problems as well, prominently by [11] for complementarity
constraints. In contrast, proximal algorithms have been developed for the
fully nonconvex structured setting, with globalization based on line search
[14, 6] and trust region [1, 3], but seem to rely on local quadratic models to
achieve sufficient decrease, hindering their direct application to (P) in
practice. Closer in spirit to our work is the Frank-Wolfe scheme of [10] with
boundable linear minimization oracle, but we proceed without any convexity
assumptions here.

#### Contributions

In this paper we investigate optimality concepts for (P) and focus on a notion
of criticality that lies between M- and B-stationarity. We propose a numerical
method based on sequential mixed-integer linearization, with a globalization
mechanism that brings together trust region and line search. The algorithm is
characterized and global convergence results are established under mild
assumptions. Finally, we report numerical results on two nontrivial example
problems, illustrating the algorithmic behavior.

##  2 Approach and Optimality Concepts

###  2.1 From projection to trust region

Attracted by the success of proximal-gradient methods to numerically solve
nonconvex nonsmooth problems, one may be tempted to address (P) with such
methods recast with sharp norms. The classical _round_ setting involves an
$\ell_{2}$-norm minimization to compute a projection onto the feasible set
$\mathcal{X}$, whereas _sharp_ projections based on $\ell_{1}$\- or
$\ell_{\infty}$-norm can be reformulated in terms of MILP. A round projected
gradient update at $\bar{x}\in\mathcal{X}$ is given by

$x^{+}\in\operatorname{proj}_{\mathcal{X}}(\bar{x}-\gamma\nabla
f(\bar{x}))=\operatorname*{\arg\min}_{x\in\mathcal{X}}\frac{1}{2\gamma}\|\bar{x}-\gamma\nabla
f(\bar{x})-x\|_{2}^{2}$

for some stepsize $\gamma>0$. This subproblem can be rewritten in view of a
quadratic model of $f$ around $\bar{x}$ as

$x^{+}\in\operatorname*{\arg\min}_{x\in\mathcal{X}}\left\\{f(\bar{x})+\langle\nabla
f(\bar{x}),x-\bar{x}\rangle+\frac{1}{2\gamma}\|x-\bar{x}\|_{2}^{2}\right\\},$
(1)

whose first-order necessary optimality conditions read

$0\in\nabla
f(\bar{x})+\frac{1}{\gamma}(x^{+}-\bar{x})+\mathcal{N}^{\textup{lim}}_{\mathcal{X}}(x^{+}),$
(2)

denoting by $\mathcal{N}^{\textup{lim}}_{\mathcal{X}}$ the limiting normal
cone to $\mathcal{X}$. Owing to the structure of $\mathcal{X}$, the round
projection requires solving a mixed-integer _quadratic_ program (MIQP),
because of the quadratic term. Seeking a MILP subproblem, one could replace
the Euclidean norm with a sharp one as distance metric and update according to

$x^{+}\in\operatorname*{\arg\min}_{x\in\mathcal{X}}\|\bar{x}-\gamma\nabla
f(\bar{x})-x\|_{p}.$ (3)

Unfortunately, taking such projections with $p\in\\{1,\infty\\}$ can lead to
odd situations. As illustrated in Figure 1 with the $\ell_{\infty}$-projection
$x_{\infty}$, not only sharp projections can be set-valued even for convex
sets, but a (badly selected) sharply projected gradient may yield an increase
of the objective, regardless of the stepsize.

Moving away from the projection point-of-view, we can consider a linear model
of $f$ around $\bar{x}$ and, to avoid artificial unboundedness, replace the
quadratic regularization term in (1) with a trust region constraint. Such
trust region subproblem at $\bar{x}\in\mathcal{X}$ with radius $\Delta>0$
reads

$\displaystyle\operatorname*{minimize}_{x\in\mathcal{X}}\quad
f(\bar{x})+\langle\nabla
f(\bar{x}),x-\bar{x}\rangle\qquad\operatorname{subject~{}to}\quad\|x-\bar{x}\|_{\textup{LP}}\leq\Delta,$
(4)

where $\|\cdot\|_{\textup{LP}}$ denotes a sharp norm on the real-valued
entries only, defined by
$\|x\|_{\textup{LP}}:=\sum_{i\notin\mathcal{I}}|x_{i}|$ and
$\|x\|_{\textup{LP}}:=\max_{i\notin\mathcal{I}}|x_{i}|$, respectively for the
$\ell_{1}$\- and $\ell_{\infty}$-norm. Despite the partial trust region
stipulation, (4) is well posed owing to Assumption (iii).

$\Delta_{1}$$\Delta_{2}$$\mathcal{X}$$\bar{x}$$-\gamma\nabla
f(\bar{x})$$x_{\infty}$$x_{2}$$x_{1}$$x_{\Delta_{1}}$$x_{\Delta_{2}}$ Figure
1: Illustration of projected gradient and trust region linearization steps
with sharp norms. A trust region step need not be single-valued nor a
projection, but satisfies the trust region constraint and yields descent.

Since the construction of (4) is similar in spirit to (1), we could expect a
connection between the two approaches, despite the different behavior
illustrated in Figure 1. For an $\ell_{p}$-projected gradient scheme, $1\leq
p\leq\infty$, the next iterate from a feasible point $\bar{x}\in\mathcal{X}$
is obtained from (3). Since $\bar{x}\in\mathcal{X}$, any solution
$x^{+}\in\mathcal{X}$ must satisfy

$\|x^{+}-\bar{x}+\gamma\nabla f(\bar{x})\|_{p}\leq\gamma\|\nabla
f(\bar{x})\|_{p}$

and, using the triangle inequality,

$\displaystyle\|x^{+}-\bar{x}\|_{p}\leq{}\|x^{+}-\bar{x}+\gamma\nabla
f(\bar{x})\|_{p}+\|-\gamma\nabla f(\bar{x})\|_{p}\leq{}2\gamma\|\nabla
f(\bar{x})\|_{p}.$

This means that, for a given feasible iterate $\bar{x}\in\mathcal{X}$, the
stepsize $\gamma>0$ controls an implicit trust region centered at $\bar{x}$,
with a direct relation between $\gamma$ and the associated
($\bar{x}$-dependent) radius $\Delta$. Then, it is clear that a projected
gradient subproblem imposes a trust region-like condition, whereas a trust
region subproblem does not seek a projection but instead to minimize (a local
model of) the objective function. Ultimately, despite the well-known
correspondence of Tikhonov regularization and $\ell_{2}$ trust region [8],
this link vanishes by considering sharp norms instead, and projected-type
updates seem to become ineffective. Thus, as on a quest to avoid round terms,
we explore the trust region perspective with sharp norms.

###  2.2 Optimality concepts

According to [14, Section 3], $\bar{x}$ is _critical_ if it is
$\gamma$-critical for some $\gamma>0$, i.e.,
$\bar{x}\in\operatorname{proj}_{\mathcal{X}}(\bar{x}-\gamma\nabla
f(\bar{x}))$, and such notion is stronger that stationarity for nonconvex
problems. Due to the difficulties discussed above in relation to round
projections, we shall consider a notion based on (4) instead. Thus, we define
the criticality measure $\Psi$ associated to (P) as follows, for all
$x\in\mathcal{X}$ and $\Delta>0$:

$\Psi(x;\Delta):=\max_{\begin{subarray}{c}w\in\mathcal{X}\\\
\|w-x\|_{\textup{LP}}\leq\Delta\end{subarray}}\langle\nabla f(x),x-w\rangle.$
(5)

Since $x\in\mathcal{X}$, the lower bound $\Psi(\cdot;\Delta)\geq 0$ holds for
all $\Delta>0$. Although analogous to the criticality measure of [16, 4],
which focused on (P) with convex $\mathcal{X}$, $\Psi$ in (5) depends
additionally on the trust region radius $\Delta$, which provides a
localization—dispensable in the convex setting. Considering (global)
minimizers of the program obtained by linearization of $f$ about $\bar{x}$ and
localization with a trust region, (5) leads to the following definition.

######  Definition 2.1 (criticality).

Given some $\Delta>0$, a point $\bar{x}\in\mathcal{X}$ is called _$\Delta$
-critical_ for (P) if $\Psi(\bar{x};\Delta)=0$; it is called _critical_ for
(P) if it is $\Delta$-critical for some $\Delta>0$.

Remark that Definition 2.1 is well-posed: if there exists some
$\bar{\Delta}>0$ such that $\Psi(\bar{x};\bar{\Delta})=0$, namely $x=\bar{x}$
is optimal with $\bar{\Delta}$, then $\bar{x}$ is necessarily optimal for all
$\Delta\in[0,\bar{\Delta}]$, for otherwise it could not be optimal for the
looser constraint with $\bar{\Delta}$. Definition 2.1 is closely related to
stationarity in nonlinear optimization, interpreted as the absence of feasible
first-order descent directions, but localized as criticality in (nonconvex)
structured optimization [14]. The B-stationarity notion of [11, Def. 2.1]
leaves out the trust region constraint and could be seen as
“$\Delta$-criticality for all $\Delta>0$”. Indeed, the simple example in
Figure 1 rules out the equivalence of B-stationarity and criticality, with the
former being stronger in general. In contrast, criticality implies
M-stationarity, based on the limiting normal cone
$\mathcal{N}^{\textup{lim}}_{\mathcal{X}}$ as in (2), whereas the converse
implication does not hold true in general. To see this fact, it suffices to
redraw Figure 1 with the negative gradient $-\nabla f(\bar{x})$ pointing
upward: while not critical, the point $\bar{x}$ can be M-stationary, by
$0\in\nabla f(\bar{x})+\mathcal{N}^{\textup{lim}}_{\mathcal{X}}(\bar{x})$.

The following Lemma 2.2 establishes that coupling constraints in $\mathcal{X}$
are essential for the hierarchy of optimality notions “B-stationarity
$\Rightarrow$ criticality $\Rightarrow$ M-stationarity”. Without such
constraints the problem in (5) can be decomposed, and the gap closed by the
integer-linearity assumption (ii) and the partial trust region norm
$\|\cdot\|_{\textup{LP}}$.

######  Lemma 2.2.

Suppose $\mathcal{X}$ does not involve constraints coupling real- and integer-
valued variables. Then, for any $\bar{x}\in\mathcal{X}$ and $\Delta>0$ it
holds that $\Psi(\bar{x},\Delta)\geq\min\\{\Delta,1\\}\Psi(\bar{x},1)$.

###### Proof.

Let $\bar{x}\in\mathcal{X}$ and $\Delta>0$ be arbitrary but fixed. The
assertion readily follows if $\Delta\geq 1$, since
$\Psi(\bar{x},\Delta)\geq\Psi(\bar{x},1)$ in that case, being the feasible set
of (5) larger. Let us focus on the case $\Delta<1$. Relying on the assumption
that $\mathcal{X}$ has no coupling constraints, we can express the feasible
set as the Cartesian product
$\mathcal{X}=\mathcal{X}_{u}\times\mathcal{X}_{z}$ for some convex set
$\mathcal{X}_{u}$ and discrete set $\mathcal{X}_{z}$. Then, from (5) we have

$\Psi(\bar{x},\Delta):={}\max_{\begin{subarray}{c}x\in\mathcal{X}\\\
\|x-\bar{x}\|_{\textup{LP}}\leq\Delta\end{subarray}}\langle\nabla
f(\bar{x}),\bar{x}-x\rangle={}\max_{\begin{subarray}{c}u\in\mathcal{X}_{u}\\\
\|u-\bar{u}\|_{\textup{LP}}\leq\Delta\end{subarray}}\langle\nabla_{u}f(\bar{x}),\bar{u}-u\rangle+\max_{z\in\mathcal{X}_{z}}\langle\nabla_{z}f(\bar{x}),\bar{z}-z\rangle,$

where the continuous component has convex feasible set and the integer term
does not depend on the trust region radius. Following the proof of [4, Lemma
3.1], this partial convexity relative to $u$ yields

$\displaystyle\Psi(\bar{x},\Delta)\geq{}$
$\displaystyle\Delta\max_{\begin{subarray}{c}u\in\mathcal{X}_{u}\\\
\|u-\bar{u}\|_{\textup{LP}}\leq
1\end{subarray}}\langle\nabla_{u}f(\bar{x}),\bar{u}-u\rangle+\max_{z\in\mathcal{X}_{z}}\langle\nabla_{z}f(\bar{x}),\bar{z}-z\rangle$
$\displaystyle\geq{}$
$\displaystyle\Delta\max_{\begin{subarray}{c}u\in\mathcal{X}_{u}\\\
\|u-\bar{u}\|_{\textup{LP}}\leq
1\end{subarray}}\langle\nabla_{u}f(\bar{x}),\bar{u}-u\rangle+\Delta\max_{z\in\mathcal{X}_{z}}\langle\nabla_{z}f(\bar{x}),\bar{z}-z\rangle$
$\displaystyle={}$ $\displaystyle\Delta\Psi(\bar{x},1)\geq{}0,$

where the second inequality is due to $\Delta<1$ and
$\bar{z}\in\mathcal{X}_{z}$. The last steps are obtained by reversing the
arguments above and using the property $\Psi(\cdot,1)\geq 0$, concluding the
proof. ∎

In the nonconvex setting considered here, the closely related Lemma 3.1 of [4]
does not hold, since $\bar{x}\in\mathcal{X}$ could be $\Delta$-critical for
some $\Delta\in(0,1)$ but, at the same time, not $1$-critical; cf. Figure 1.
Thus, if possible at all, an analogous result including coupling constraints
must be restricted to the noncritical case, as it cannot hold in general for
critical points. Furthermore, the stipulation of a _connected_ feasible set
$\mathcal{X}$ does not resolve the gap between criticality and B-stationarity.
Comparing with [11], we conjecture that the complementarity structure there
allowed to close the gap and dispose of the localization in (5), but this is
not possible in general.

Before developing and analysing a numerical scheme for (P), we shall introduce
an approximate counterpart of Definition 2.1. The following definition
requires that a point $\bar{x}\in\mathcal{X}$ is not too far from being a
minimizer for (4) in terms of function value, as monitored using $\Psi$ from
(5).

######  Definition 2.3 ($\varepsilon$-criticality).

Given some $\varepsilon>0$ and $\Delta>0$, a point $\bar{x}\in\mathcal{X}$ is
called _$\varepsilon$ -$\Delta$-critical_ for (P) if
$\Psi(\bar{x};\Delta)\leq\varepsilon$; it is called _$\varepsilon$ -critical_
for (P) if it is $\varepsilon$-$\Delta$-critical for some $\Delta>0$.

For later reference we give the following characterization of noncritical
points.

######  Lemma 2.4.

Let $\bar{x}\in\mathcal{X}$ be not critical for (P). Then, $\bar{x}$ is not
isolated, $\nabla f(\bar{x})\neq 0$ and, for all $\Delta>0$, there exist
$x\in\mathcal{X}$ and $\varepsilon>0$ such that
$\|x-\bar{x}\|_{\textup{LP}}\leq\Delta$ and $\langle\nabla
f(\bar{x}),x-\bar{x}\rangle\leq-\varepsilon$.

###### Proof.

If $\bar{x}$ is isolated or $\nabla f(\bar{x})=0$, then
$\bar{x}\in\mathcal{X}$ must be critical according to Definition 2.1, hence
the first two assertions. The third is obtained similarly, by negating the
criticality property for $\bar{x}$; see (5). ∎

##  3 Sequential mixed-integer linearization method

We now introduce our iterative numerical scheme for (P) based on mixed-integer
linear approximations. Then, we will investigate its behavior and show that it
converges to critical points under mild assumptions.

###  3.1 Algorithm

The proposed method, detailed in Algorithm 3.1, builds upon a sequence of
subproblems of the form (4). Since MILP subproblems take most of the
computational toll, with the remaining steps requiring only negligible effort,
the problem structure as well as modern MILP technology should be exploited to
limit the computational burden; see [11, 10]. After the trust region
linearized subproblem at Algorithm 3.1, which guarantees that
$x^{k+1}\in\mathcal{X}$, $\|x^{k+1}-x^{k}\|_{\textup{LP}}\leq\Delta_{k}$, and

$\Psi(x^{k};\Delta_{k})=\Psi_{k}:=\langle\nabla
f(x^{k}),x^{k}-x^{k+1}\rangle\geq 0,$ (6)

the quality of the tentative update $x^{k+1}$ from $x^{k}$ is tested. In
particular, we monitor criticality and terminate at Algorithm 3.1 if
(approximately) satisfied. Then, a sufficient decrease condition requires that

$f(x^{k+1})\leq\Phi_{k}-\varrho\Psi(x^{k};\Delta_{k})$ (7)

holds for the update $x^{k+1}$ to be accepted; otherwise, similarly to a
backtracking line search procedure, the trust region radius $\Delta_{k}$ is
reduced and the attempt discarded. With the lens of trust region methods, (7)
compares the predicted linearized improvement againt the actual change, with
the ratio parameter $\varrho>0$ acting as threshold for acceptance. In fact,
the actual change at Algorithm 3.1 refers to a merit value $\Phi_{k}$ which
coincides with the cost only if $p=1$, owing to Algorithm 3.1. For $p\in(0,1)$
the acceptance criterion is less conservative and potentially allows larger
radii and faster convergence in practice, with the objective $f(x^{k})$
decreasing nonmonotonically; see [5].

Input: $x^{0}\in\mathcal{X}$, $\varepsilon>0$

Data: $\Delta_{0}>0$, $\varrho,\kappa\in(0,1)$, $p\in(0,1]$

1 set $\Phi_{0}\leftarrow f(x^{0})$

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

3 $x^{k+1}\in\operatorname*{\arg\min}_{x\in\mathcal{X}}\langle\nabla
f(x^{k}),x\rangle\colon~{}\|x-x^{k}\|_{\textup{LP}}\leq\Delta_{k}$

4 $a_{k}\leftarrow\Phi_{k}-f(x^{k+1})$

5 $\Psi_{k}\leftarrow\langle\nabla f(x^{k}),x^{k}-x^{k+1}\rangle$

6 if  _$\Psi_{k}\leq\varepsilon$_ then return $x^{k}$ //
$\varepsilon$-$\Delta_{k}$-critical

7 if  _$a_{k} <\varrho\Psi_{k}$_ then

8 $\Delta_{k}\leftarrow\kappa\Delta_{k}$ and go to Algorithm 3.1

9 set $\Phi_{k+1}\leftarrow(1-p)\Phi_{k}+pf(x^{k+1})$

10 select $\Delta_{k+1}$ by trust region update (8) or reset (9)

11

Algorithm 3.1 Successive mixed-integer linearization algorithm for (P)

Finally, the radius is updated at Algorithm 3.1, where one can apply classical
trust region rules such as

$\Delta_{k+1}:=\begin{cases}\kappa\Delta_{k}&\text{if}~{}\rho_{k}<\varrho_{1},\\\
\Delta_{k}&\text{if}~{}\varrho_{1}\leq\rho_{k}<\varrho_{2},\\\
\Delta_{k}/\kappa&\text{if}~{}\varrho_{2}\leq\rho_{k}\\\ \end{cases}$ (8)

where $\rho_{k}:=a_{k}/\Psi_{k}$, or the line search-like reset

$\Delta_{k+1}\in[\Delta_{\min},\Delta_{\max}],$ (9)

given parameters such that $0<\varrho\leq\varrho_{1}<\varrho_{2}<1$ or
$0<\Delta_{\min}\leq\Delta_{\max}$. Remark that considering iteration-
dependent algorithmic parameters would be possible too.

We note that the input requirement in Algorithm 3.1 of a feasible initial
point is not a restriction: given some $x^{0}\notin\mathcal{X}$, it is enough
to (sharply) project it onto the feasible set by solving a MILP, for instance
$\min_{x\in\mathcal{X}}\|x-x^{0}\|_{1}$, and run the algorithm from a solution
thereof.

An optional refinement step can also be integrated in Algorithm 3.1, similarly
to EQP steps of SLQP approaches [4, 8] and to BQP steps of [11]. Here we can
consider NLP steps, by fixing the integer-valued variables and solving a
linearly-constrained nonlinear program (NLP). Without impairing the
convergence properties, this optional step can improve the objective value
while maintaining feasibility, taking advantage of efficient NLP methods and
possibly exploiting additional smoothness of $f$, as discussed in Section 4
below.

Finally, the termination criterion at Algorithm 3.1 is based on the
approximate criticality notion of Definition 2.3, and supported by the
following convergence analysis, which guarantees finite termination for any
$\varepsilon>0$.

###  3.2 Convergence analysis

To derive our convergence results, we further assume that

  1. (iv)

function $f$ is continuously differentiable and $\nabla f$ is locally
Lipschitz continuous.

Equipped with a preliminary result, we show in Lemma 3.2 that the inner loop
of Algorithm 3.1 always terminates after finitely many iterations.
Consequently, Algorithm 3.1 is well-defined and produces a sequence of
feasible iterates $\\{x^{k}\\}$ along with decreasing merit values; cf. Lemma
3.3. Moreover, all accumulation points of the sequence $\\{x^{k}\\}$ are
feasible for (P) because $\mathcal{X}$ is closed and $x^{k}\in\mathcal{X}$ for
each $k\in\mathbb{N}$. Then, it remains to prove that the accumulation points
are also critical. After the asymptotic characterization in Lemma 3.4, Lemma
3.5 shows that, in a neighborhood of a feasible but noncritical point, the
subproblem at Algorithm 3.1 eventually generates an acceptable step while the
trust region radius remains bounded away from zero. Finally, Theorem 3.6
synthesizes the main convergence results.

We start with a preliminary result about the reduction implied by Algorithm
3.1, which follows from Taylor’s theorem. Notice how Assumption (ii) allows to
drop the integer-valued variables from the quadratic term of the expansion.

######  Lemma 3.1.

Consider a function $f$ under the standing assumptions and a nonempty closed
bounded set $\Omega\subseteq\operatorname{dom}f$. Partitioning $x=(u,z)$ with
$u$ and $z$ the real- and integer-valued variables, respectively, define

$L:=\sup\left\\{\frac{\|\nabla f(x_{1})-\nabla
f(x_{2})\|_{1}}{\|u_{1}-u_{2}\|_{\infty}}\,\middle|\,\begin{aligned}
&x_{1},x_{2}\in\Omega,\\\ &u_{1}\neq u_{2}\end{aligned}\right\\}.$

Then $L<\infty$, and for all $x_{1},x_{2}\in\Omega$ the upper bound

$f(x_{2})\leq f(x_{1})+\langle\nabla
f(x_{1}),x_{2}-x_{1}\rangle+\frac{L}{2}\|u_{2}-u_{1}\|_{\infty}^{2}$

is valid.

###### Proof.

By Assumption (ii) it is $\nabla f(x)=(\nabla f_{1}(u),f_{2})$, hence
$L<\infty$ coincides with the Lipschitz constant of $\nabla f$ over $\Omega$,
finite by Assumption (iv). Then, the asserted inequality follows from Taylor-
type arguments, by considering the variation from $x_{1}$ to $x_{2}$ along
real- and integer-valued components separately.

We sketch the second part of the proof as follows. Let us define
$\widehat{x}:=(u_{2},z_{1})$ and, for simplicity, assume
$\widehat{x}\in\Omega$. Using Assumption (ii) and then Taylor’s theorem, e.g.
[11, Lemma 4.3], we obtain the inequality

$f(x_{2})=f(\widehat{x})+\langle\nabla
f(\widehat{x}),x_{2}-\widehat{x}\rangle\leq f(x_{1})+\langle\nabla
f(x_{1}),\widehat{x}-x_{1}\rangle+\frac{L}{2}\|\widehat{x}-x_{1}\|_{\infty}^{2}+\langle\nabla
f(\widehat{x}),x_{2}-\widehat{x}\rangle.$

Finally, exploiting the structure of $\nabla f$ and the definition of
$\widehat{x}$, one recognizes the asserted upper bound. ∎

Next, we employ Lemma 3.1 to prove that the iterations of Algorithm 3.1 are
well-defined.

######  Lemma 3.2 (well-definedness).

Consider (P) under the standing assumptions. Then the inner loop of Algorithm
3.1 terminates in a finite number of steps.

###### Proof.

If $x^{k}$ is critical, then it is $\bar{\Delta}$-critical for some
$\bar{\Delta}>0$, and $x=x^{k}$ solves the subproblem at Algorithm 3.1 as soon
as $\Delta_{k}\leq\bar{\Delta}$. Regardless of the choice
$x^{k+1}\in\mathcal{X}$, it must yield $\Psi_{k}=0$ and immediate termination
of the algorithm, for any $\varepsilon>0$.

If $x^{k}$ is not critical, then it must be that $\nabla f(x^{k})\neq 0$ by
Lemma 2.4. The claim will follow by induction since $\Phi_{0}=f(x^{0})$ and,
for each $k\in\mathbb{N}$, the acceptance according to Algorithm 3.1 and the
update rule at Algorithm 3.1 imply that

$\Phi_{k+1}:={}(1-p)\Phi_{k}+pf(x^{k+1})\geq{}(1-p)\left[f(x^{k+1})+\varrho\Psi_{k}\right]+pf(x^{k+1})\geq{}f(x^{k+1}),$
(10)

proving that the nonmonotone acceptance criterion is less conservative than
the monotone one (which has $p=1$ and $\Phi_{k}=f(x^{k})$ for all
$k\in\mathbb{N}$). Moreover, all attempts satisfy $x^{k+1}\in\mathcal{X}$ and
$\|x^{k+1}-x^{k}\|_{\textup{LP}}\leq\Delta_{k}$. In particular, all attempts
$x^{k+1}$ remain in a bounded set, owing to Assumption (iii), so there exists
a constant $L_{k}>0$ related to $\nabla f$ over that set, in the sense of
Lemma 3.1. Then, all attempts satisfy the upper bound

$f(x^{k+1})\leq f(x^{k})+\langle\nabla
f(x^{k}),x^{k+1}-x^{k}\rangle+\frac{L_{k}}{2}\|u^{k+1}-u^{k}\|_{\infty}^{2}.$

Combined with (10) and then (6), this yields

$\displaystyle\rho_{k}:=\frac{a_{k}}{\Psi_{k}}\geq{}\frac{f(x^{k})-f(x^{k+1})}{\Psi_{k}}\geq{}\frac{\Psi_{k}-\frac{L_{k}}{2}\|u^{k+1}-u^{k}\|_{\infty}^{2}}{\Psi_{k}}\geq{}1-\frac{L_{k}}{2}\frac{\|u^{k+1}-u^{k}\|_{\infty}^{2}}{\|\nabla
f(x^{k})\|\|x^{k+1}-x^{k}\|}$

where the last step follows from the triangle inequality. Due to the
equivalence of norms in $\mathbb{R}^{n}$ there exists a constant
$c_{\infty}>0$ such that for each $d\in\mathbb{R}^{n}$ it holds that
$c_{\infty}\|d\|_{\infty}\leq\|d\|$. Similarly, there is a constant
$c_{\textup{LP}}>0$ relating $\|u\|_{\infty}$ and $\|x\|_{\textup{LP}}$.
Therefore, since $\|x^{k+1}-x^{k}\|\geq\|u^{k+1}-u^{k}\|$, we have that

$\displaystyle\rho_{k}\geq{}1-\frac{L_{k}}{2}\frac{\|u^{k+1}-u^{k}\|_{\infty}^{2}}{\|\nabla
f(x^{k})\|\|u^{k+1}-u^{k}\|}\geq{}1-\frac{L_{k}}{2c_{\infty}}\frac{\|u^{k+1}-u^{k}\|_{\infty}}{\|\nabla
f(x^{k})\|}\geq{}1-\frac{L_{k}c_{\textup{LP}}}{2c_{\infty}}\frac{\Delta_{k}}{\|\nabla
f(x^{k})\|}$

where the last step uses the trust region stipulation with radius
$\Delta_{k}$. Since $L_{k}$ and $x^{k}$ are fixed, with $\|\nabla
f(x^{k})\|>0$, while $\Delta_{k}$ is decreased, the condition
$\rho_{k}\geq\varrho$ becomes satisfied as soon as

$\Delta_{k}\leq\frac{2c_{\infty}}{L_{k}c_{\textup{LP}}}(1-\varrho)\|\nabla
f(x^{k})\|$ (11)

and the inner loop stops at Algorithm 3.1, after finitely many backtracks
since $\kappa\in(0,1)$. ∎

######  Lemma 3.3 (descent behavior).

Consider (P) under the standing assumptions and the iterates generated by
Algorithm 3.1. Then,

  1. (i)

the sequence $\\{\Phi_{k}\\}$ is monotonically decreasing and, for every
$k\in\mathbb{N}$, one has $\Psi_{k}\geq 0$ and

$f(x^{k+1})+(1-p)\varrho\Psi_{k}\leq\Phi_{k+1}\leq\Phi_{k}-p\varrho\Psi_{k}.$
(12)

  2. (ii)

Every iterate $x^{k}$ is feasible and remains in the sublevel set
$\\{x\in\mathcal{X}\,|\,f(x)\leq f(x^{0})\\}$.

###### Proof.

Based on Lemma 3.2, the iterates of Algorithm 3.1 are well-defined. For (i),
by the update rule at Algorithm 3.1 and the condition at Algorithm 3.1, we
have that

$\displaystyle\Phi_{k+1}:={}(1-p)\Phi_{k}+pf(x^{k+1})\leq{}(1-p)\Phi_{k}+p\left[\Phi_{k}-\varrho\Psi_{k}\right]={}\Phi_{k}-p\varrho\Psi_{k}$

holds upon acceptance. The lower bound on $\Phi_{k}$ is readily obtained from
(10). That $\Psi_{k}\geq 0$ for each $k\in\mathbb{N}$ follows from (6).

Consider (ii). Owing to $x^{0}\in\mathcal{X}$ and continuity of $f$, the
sublevel set of $f$ induced by $x^{0}$ is a subset of $\mathcal{X}$. Then,
since $p\in(0,1]$ and $\varrho\Psi_{k}\geq 0$, it follows from
$\Phi_{0}:=f(x^{0})$ and (i) that
$f(x^{k})\leq\Phi_{k}\leq\ldots\leq\Phi_{0}:=f(x^{0})<\infty$, whence the
inclusion of $x^{k}$ in the sublevel set. ∎

######  Lemma 3.4.

Consider (P) under the standing assumptions and suppose
$\inf_{\mathcal{X}}f>-\infty$. Then, relative to the iterates generated by
Algorithm 3.1,

  1. (i)

$\\{f(x^{k})\\}$ and $\\{\Phi_{k}\\}$ converge, the latter from above, to a
finite value $f_{\star}\geq\inf_{\mathcal{X}}f$;

  2. (ii)

$\sum_{k\in\mathbb{N}}\Psi_{k}<\infty$;

  3. (iii)

$\lim_{k\to\infty}\Psi_{k}=0$.

###### Proof.

Based on Lemma 3.2, the iterates of Algorithm 3.1 are well-defined and we may
assume the sequence $\\{x^{k}\\}$ is infinite (it suffices to skip Algorithm
3.1). The assertion in (i) regarding $\\{\Phi_{k}\\}$ follows from (12) with
the premises. Then, considering the update rule at Algorithm 3.1,

$pf(x^{k+1})=\Phi_{k+1}-(1-p)\Phi_{k}=p\Phi_{k}+(\Phi_{k+1}-\Phi_{k}),$

the convergence of $\\{\Phi_{k}\\}$ implies that of $\\{f(x^{k})\\}$ to the
same value $f_{\star}$, since $p\in(0,1]$.

Regarding (ii), a telescoping argument on (12), together with (i), yields

$\Phi_{0}-f_{\star}\geq\sum_{k=0}^{j}\Phi_{k}-\Phi_{k+1}\geq
p\sum_{k=0}^{j}\delta_{k+1}=p\varrho\sum_{k=0}^{j}\Psi_{k}\geq 0.$

The finite sum follows from the independence of the (finite) upper bound
$\Phi_{0}-f_{\star}$ on $j\in\mathbb{N}$.

Finally, assertion (iii) follows from (ii), since $\Psi_{k}\geq 0$ by Lemma
3.3. ∎

Owing to (6), Lemma 3.4(iii) gives that $\\{\Psi(x^{k};\Delta_{k})\\}\to 0$ if
$\inf_{\mathcal{X}}f>-\infty$. Therefore, if $\\{\Delta_{k}\\}$ remains
bounded away from zero, we expect accumulation points of $\\{x^{k}\\}$ to be
critical. Indeed, for any noncritical point $x^{\infty}\in\mathcal{X}$ and
sequence $\\{x^{k}\\}$ such that $x^{k}\to x^{\infty}$, the trust region
radius $\\{\Delta_{k}\\}$ can be bounded from below by some positive constant.
In the following result, inspired by [11, Lemma 4.5], by _relative
neighborhood_ we mean the intersection of a neighborhood of $x^{\infty}$ in
$\mathbb{R}^{n}$ with the feasible set $\mathcal{X}$.

######  Lemma 3.5.

Consider (P) under the standing assumptions. Let a point
$x^{\infty}\in\mathcal{X}$ be not critical for (P). Then, there exists a
relative neighborhood $N^{\infty}$ of $x^{\infty}$ such that for any sequence
$\\{x^{k}\\}\subset N^{\infty}$ with $x^{k}\to x^{\infty}$ the steps of
Algorithm 3.1 generate a sequence $\\{\Delta_{k}\\}$ bounded away from zero.

###### Proof.

Because $x^{\infty}$ is not critical, Lemma 2.4 ensures the existence of a
nonempty non-singleton relative neighborhood $N^{\infty}$ of $x^{\infty}$ and
that $\nabla f(x^{\infty})\neq 0$. Now, by continuity of $\nabla f$ we can
find $N^{\infty}$ such that $\inf_{N^{\infty}}\|\nabla f\|\geq\nabla
f_{\infty}$ for some $\nabla f_{\infty}>0$. Denote $L_{\infty}$ the Lipschitz
constant of $\nabla f$ over $N^{\infty}$, finite by Assumption (iv). Since
$\Delta_{k}$ is reduced when needed for sufficient decrease but not otherwise,
$\Delta_{k}$ cannot get lower than (a fixed constant or a fraction of) the
worst-case threshold for acceptance at Algorithm 3.1, by Lemma 3.2. Precisely,
from (11) we have the lower bound
$\Delta_{k}\geq\min\\{\Delta_{\min},\kappa\Delta_{\infty}\\}$, where
$\Delta_{\infty}:=2(1-\varrho)c_{\infty}\nabla
f_{\infty}/(L_{\infty}c_{\textup{LP}})>0$. ∎

By virtue of Lemma 3.5, the trust region radius does not vanish in the
vicinity of noncritical points, hence sufficiently small steps can always be
accepted there, which deliver sufficient reduction and improvement. We are now
in a position to prove our main convergence result, connecting accumulation
points and criticality. The following Theorem 3.6 is akin to [4, Thm 3.8] and
[11, Thm 4.2].

######  Theorem 3.6.

Consider (P) under the standing assumptions and let a point
$x^{0}\in\mathcal{X}$ be arbitrary but fixed. Then, relative to the iterates
generated by Algorithm 3.1 with $\varepsilon=0$, one of the following mutually
exclusive outcomes must occur:

  1. (i)

$\lim_{k\to\infty}f(x^{k})=-\infty$;

  2. (ii)

Algorithm 3.1 terminates at a critical point, that is, $x^{k}$ solves the
subproblem at Algorithm 3.1 for some $k\in\mathbb{N}$;

  3. (iii)

Algorithm 3.1 generates an infinite sequence of iterates $\\{x^{k}\\}$ with
decreasing merit values $\\{\Phi_{k}\\}$. If $\\{x^{k}\\}$ has an accumulation
point $x^{\infty}$, then $x^{\infty}$ is feasible and critical for (P).

###### Proof.

We need to consider only the last outcome because in the other cases we either
observe lower unboundedness or obtain a critical point by virtue of the
termination condition.

Since the inner loop terminates finitely by Lemma 3.2, and each attempt is
feasible for the subproblem at Algorithm 3.1, we obtain by induction over the
iterations that $x^{k}\in\mathcal{X}$ for all $k\in\mathbb{N}$. Furthermore,
by Lemma 3.3 all iterations yield an improvement in terms of merit —and the
objective value eventually decreases as well, owing to Lemma 3.4. It remains
to show that every accumulation point of $\\{x^{k}\\}$ is critical.

Seeking a contradiction, we assume that an accumulation point $x^{\infty}$ of
the sequence $\\{x^{k}\\}$ is not critical. Let us denote $\\{x^{k}\\}_{k\in
K}$ a subsequence such that $x^{k}\to_{K}x^{\infty}$ and, possibly relabeling,
restrict $\\{x^{k}\\}_{k\in K}$ to the neighborhood $N^{\infty}$ defined in
Lemma 3.5, so that $x^{k}\in N^{\infty}$ for all $k\in K$.

Now we show that $\\{\Psi_{k}\\}_{k\in K}$ remains bounded away from zero.
First of all, by noncriticality of $x^{\infty}$, Lemma 2.4 gives that
$\langle\nabla f(x^{\infty}),x^{\infty}-x\rangle\geq\varepsilon_{0}$ for some
$x\in\mathcal{X}$ and $\varepsilon_{0}>0$. Then, recall the continuity of
$\nabla f$, the minimality property of $x^{k+1}$ for $\langle\nabla
f(x^{k}),\cdot-x^{k}\rangle$ within a radius $\Delta_{k}$ (bounded away from
zero by Lemma 3.5), and the convergence $x^{k}\to x^{\infty}$. These facts
imply the existence of some $\varepsilon>0$ such that
$\Psi_{k}\geq\varepsilon$ holds for all $k\in K$.

Combined with (12), and for any $k_{0}\in K$, the boundedness of
$\\{\Psi_{k}\\}_{k\in K}$ away from zero results in

$f(x^{k})\leq{}\Phi_{k}={}\Phi_{k_{0}}+\sum_{j=k_{0}}^{k-1}(\Phi_{j+1}-\Phi_{j})\leq{}\Phi_{k_{0}}-\sum_{j=k_{0}}^{k-1}p\varrho\Psi_{j}\leq{}\Phi_{k_{0}}-\sum_{j=k_{0}}^{k-1}p\varrho\varepsilon={}\Phi_{k_{0}}-(k-k_{0})p\varrho\varepsilon\to-\infty$

as $k\to_{K}\infty$. This fact contradicts the convergence
$f(x^{k})\to_{K}f(x^{\infty})$ by $x^{k}\to_{K}x^{\infty}$, thus showing that
every accumulation point $x^{\infty}$ of $\\{x^{k}\\}$ is critical. ∎

Finally, in analogy with Theorem 3.6, finite termination with an
$\varepsilon$-critical point can be established for any $\varepsilon>0$, if
$\inf_{\mathcal{X}}f>-\infty$, via Lemma 3.4.

##  4 Numerical Examples

We report on numerical results for two example problems: the discrete-time
optimal control of a hybrid system with hysteresis and the concurrent design
and operation planning of a processing network.

To exploit the problem structure and the available technology, we consider
also a variant of Algorithm 3.1 with the NLP refinement step suggested in
Section 3.1. In practice, when the integer-valued variables remain untouched
by a successful iteration (namely, after sufficient decrease is attained), the
real-valued variables are polished by solving the nonlinear program obtained
from (P) by fixing the integer-valued variables at their current value. More
precisely, the NLP refinement may take place before executing Algorithm 3.1,
only if $z^{k+1}=z^{k}$. An improved $u^{k+1}$ is obtained by solving (up to
$\varepsilon$-stationarity) the linearly-constrained minimization problem
arising from (P) with the additional constraint $z=z^{k}$.

The numerical results are generated with a prototype MATLAB implementation of
Algorithm 3.1, adopting intlinprog for the MILPs arising at Algorithm 3.1 and
fmincon for the optional NLPs. The default algorithmic parameters are set as
follows: tolerance $\varepsilon=10^{-8}$, initial radius $\Delta_{0}=1$,
sufficient decrease parameter $\varrho=\nicefrac{{1}}{{10}}$, reduction factor
$\kappa=\nicefrac{{1}}{{2}}$, and monotonicity parameter
$p=\nicefrac{{1}}{{2}}$. Algorithm 3.1 uses the trust region update rule (8)
with $\varrho_{1}=\varrho$ and $\varrho_{2}=2\varrho$. If an infeasible
starting point $x^{0}$ is provided, a feasible one is computed by solving
$\min_{x\in\mathcal{X}}\|x-x^{0}\|_{1}$ reformulated as a MILP. Moreover, the
execution of Algorithm 3.1 is stopped if the MILP solver fails, either by
declaring the subproblem at Algorithm 3.1 infeasible or resulting in a
negative criticality measure $\Psi_{k}$ at Algorithm 3.1. These occurrences
are both due to numerical issues within the MILP solver, as the subproblem is
guaranteed to be feasible and the criticality measure to be nonnegative, by
construction. All source codes used to generate the results in Section 4 have
been archived on Zenodo at doi: 10.5281/zenodo.10007958.

###  4.1 Car with hysteretic turbo charger

In this section we apply Algorithm 3.1 to a numerical example for the point-
to-point optimal control problem of a car with turbo charger. Inspired by [12,
Section V], we consider a double-integrator point mass model equipped with a
turbo accelerator subject to hysteresis effects.

The car is described by its position $q(t)$, velocity $v(t)$ and turbo state
$w(t)\in\\{0,1\\}$ for each time $t\in[0,T]$, $T>0$. The control variables are
the acceleration $a(t)$ and brake $b(t)$ pedals. The turbo is activated when
the velocity exceeds $v\geq v_{+}:=10$ and is deactivated when it falls below
$v\leq v_{-}:=5$; when it is on, it makes the nominal thrust three times more
effective. In contrast, the braking force remains unaffected by the turbo
state.

In summary, the state vector reads as $(q,v,w)\in\mathbb{R}^{3}$ with dynamics
$\dot{q}=v$, $\dot{v}=f-b$, where the thrust $f$ has two modes of operation,
depending on the turbo state, described by $f(a)=a$ if $w=0$ and $f(a)=3a$ if
$w=1$. The turbo state $w$ behaves according to the hysteresis characteristic
described above. The acceleration control is bounded by $0\leq a\leq
a_{\max}$, brake control by $0\leq b\leq b_{\max}$, and the velocity by
$|v|\leq v_{\max}$. The initial state is $q(0)=v(0)=w(0)=0$. Given a final
time $T>0$ and position $q_{\text{end}}>0$ and parameters
$\alpha_{a},\alpha_{b}\geq 0$, we are interested in minimizing the control
effort
$\alpha_{a}\int_{0}^{T}u^{2}(t)\mathrm{d}t+\alpha_{b}\int_{0}^{T}b^{3}(t)\mathrm{d}t$
to achieve $q(T)=q_{\text{end}}$ and $v(T)=0$.

By direct discretization [2, Section 4.2] one can formulate this optimal
control problem with hybrid dynamics as a (linear) GDP. Then, rewriting
logical propositions as big-M constraints results in the form (P). Let us
consider a time grid $0=t_{0}<\ldots<t_{N}=T$, with $N$ intervals of size
$h=T/N$, and state and control approximations there, denoted
$q_{k},v_{k},\ldots,b_{k}$ and such that $q_{k}\approx q(t_{k})$ for each
$k=0,\ldots,N$. Auxiliary variables $f_{k}$ for the thrust $f$ are included
for clarity. Then, the integral objective is approximated via the trapezoidal
rule by

$\alpha_{a}h\sum_{k=0}^{N-1}\frac{a_{k}^{2}+a_{k+1}^{2}}{2}+\alpha_{b}h\sum_{k=0}^{N-1}\frac{b_{k}^{3}+b_{k+1}^{3}}{2}$
(13)

and so are the dynamics with finite differences, for $k=0,\ldots,N-1$,

$\displaystyle\frac{q_{k+1}-q_{k}}{h}=\frac{v_{k+1}+v_{k}}{2},$
$\displaystyle\frac{v_{k+1}-v_{k}}{h}=\frac{f_{k+1}+f_{k}}{2}-\frac{b_{k+1}+b_{k}}{2}.$
(14)

Here, the thrust $f_{k}$ can be encoded by a linear big-M model, instead of
_if-else_ cases, as

$\displaystyle f_{k}-u_{k}\leq{}Mw_{k},$ $\displaystyle
f_{k}-u_{k}\geq{}-Mw_{k},$ (15) $\displaystyle f_{k}-3u_{k}\leq{}M(1-w_{k}),$
$\displaystyle f_{k}-3u_{k}\geq{}-M(1-w_{k})$

where $M>0$ is fixed and large enough. The hysteresis characteristic
corresponds to logical propositions describing activation and deactivation
conditions, in the form

$\displaystyle w_{k}=0~{}\wedge~{}v_{k}>v_{+}~{}\Rightarrow~{}w_{k+1}=1,$
$\displaystyle w_{k}=1~{}\wedge~{}v_{k}<v_{-}~{}\Rightarrow~{}w_{k+1}=0,$
$\displaystyle w_{k}=0~{}\wedge~{}v_{k}<v_{+}~{}\Rightarrow~{}w_{k+1}=0,$
$\displaystyle w_{k}=1~{}\wedge~{}v_{k}>v_{-}~{}\Rightarrow~{}w_{k+1}=1$

for $k=0,\ldots,N-1$, with $\vee$ denoting the logical “and”. Rewriting as
clauses and using De Morgan’s laws, these too can be cast as algebraic big-M
constraints, reading respectively

$\displaystyle v_{k}\leq v_{+}+M(w_{k}+w_{k+1}),$ $\displaystyle v_{k}\geq
v_{-}-M(2-w_{k}-w_{k+1}),$ (16) $\displaystyle v_{k}\geq
v_{+}-M(w_{k}+1-w_{k+1}),$ $\displaystyle v_{k}\leq v_{-}+M(1-w_{k}+w_{k+1}).$

Finally, bounds are imposed only on the time grid, analogously to boundary
conditions, as

$\displaystyle 0\leq a_{k}\leq a_{\max},$ $\displaystyle 0\leq b_{k}\leq
b_{\max},$ $\displaystyle-v_{\max}\leq v_{k}\leq v_{\max}$
$\displaystyle\text{for}~{}k=0,\ldots,N,$ (17) $\displaystyle
q_{0}=v_{0}=w_{0}=0,$ $\displaystyle q_{N}=q_{\text{end}},$ $\displaystyle
v_{N}=0.$

Overall, a discretization with $N$ intervals results in a problem (P)
according to (13)–(17) with $5(N+1)$ real-valued and $(N+1)$ binary-valued
decision variables, $2N$ linear equality constraints, $4(N+1)+4N$ linear
inequality constraints, and simple bounds.

#### Setup

Consider the problem instances with $N\in\\{25,50,100\\}$ discretization
intervals and parameters $\alpha_{a}=1$, $\alpha_{b}=10^{-2}$,
$q_{\text{end}}=150$, $T=10$, $u_{\max}=5$, $b_{\max}=10$, $v_{\max}=25$, and
$M=20$. For each instance we run Algorithm 3.1 starting from 100 initial
points obtained by sampling a normal distribution with zero mean and standard
deviation 10.

#### Results

Algorithm 3.1 successfully returned in all cases with a critical point (within
the specified tolerance). Solutions are depicted as trajectories in Figure 2
and computational results are summarized in Table 1. Out of 100 runs, at most
4 different solutions are returned by Algorithm 3.1 for each instance. These
have essentially the same pattern but slightly perturbed turbo activation
point, as shown by Figure 2 for $N=100$. Moreover, Table 1 highlights that not
only all solutions have very similar objective value, but also that the
solution process is relatively consistent. In particular, the number of
iterations and successfull updates are in narrow ranges, independent on the
discretization level. Overall, despite the straightforward implementation and
problem size ($n=606$ when $N=100$), Algorithm 3.1 was able to achieve
accurate results in few seconds.

$0$$2$$4$$6$$8$$10$$0$$50$$100$$150$$t$$q$$0$$2$$4$$6$$8$$10$$0$$10$$20$$t$$v$$0$$2$$4$$6$$8$$10$$-10$$0$$10$$t$$f-b$ Figure 2: Solutions for the problem of a car with turbo charger with $N=100$ discretization intervals: superimposed states and control trajectories of all test runs, with thresholds (black solid lines) and different markers when the turbo is active (red circle) or inactive (blue dot). |  | min | q(25%) | median | q(75%) | max
---|---|---|---|---|---|---
$N=25$ | objective value | 71.20 | 71.20 | 71.20 | 71.20 | 74.35
| iterations | 125 | 130 | 131 | 133 | 141
| successful updates | 59 | 61 | 63 | 64 | 69
| runtime [s] | 0.44 | 0.58 | 0.62 | 0.67 | 0.91
|  $\operatorname{proj}(x^{0})$ runtime [s] | 0.06 | 0.13 | 0.15 | 0.18 | 0.29
$N=50$ | objective value | 69.46 | 69.80 | 69.80 | 69.80 | 71.74
| iterations | 128 | 134 | 136 | 137 | 163
| successful updates | 60 | 65 | 66 | 68 | 84
| runtime [s] | 1.01 | 1.30 | 1.37 | 1.44 | 1.93
|  $\operatorname{proj}(x^{0})$ runtime [s] | 0.26 | 0.56 | 0.66 | 0.76 | 1.71
$N=100$ | objective value | 68.66 | 68.95 | 69.6 | 69.6 | 69.6
| iterations | 134 | 136 | 137 | 139 | 141
| successful updates | 66 | 68 | 70 | 71 | 73
| runtime [s] | 2.84 | 3.26 | 3.43 | 3.53 | 4.34
|  $\operatorname{proj}(x^{0})$ runtime [s] | 1.40 | 3.03 | 3.82 | 7.49 | 124.08
Table 1: Summary of computational results for the problem of a car with turbo
charger. Statistics (minimum, maximum, and quartiles) over different initial
points. Algorithm 3.1’s runtime does not include the sharp projection of
infeasible initial points, denoted $\operatorname{proj}(x^{0})$.

###  4.2 Processing network design and operation

In this section we consider the synthesis of a processing network that was
originally formulated as a GDP in [15, Example 2]. The problem entails the
concurrent design and operation planning of a networked system, that is, nodes
and flows along the edges are jointly optimized. Figure 3 shows the
superstructure which involves the possible selection of 8 interconnected
processing units. With this illustrative example we demonstrate the potential
benefits of including NLP refinement steps in Algorithm 3.1.

The original model consists of 25 real-valued $u$ and 8 Boolean variables $Z$
to model flows on the lines and the existence or non-existence of processing
units. The objective function includes costs for the selected processing
units, operating costs and revenues from sales of products. The formulation
involves material balance and flow constraints, logical propositions, and
disjunctions for consistent physical modelling. Boolean specifications are
transformed into linear algebraic constraints on binary variables; for
instance, the proposition $Z_{3}~{}\Rightarrow~{}Z_{1}\vee Z_{2}$ is rewritten
as the clause $Z_{1}\vee Z_{2}\vee\neg Z_{3}$ and then transformed into
$z_{1}+z_{2}\geq z_{3}$. Specifications such as $Z_{1}\veebar Z_{2}$, with
$\veebar$ the logical “xor”, are encoded as algebraic linear constraints on
the binary-valued variables: $z_{1}+z_{2}=1$; see also [15, Table 4].
Disjunctions involving linear constraints, e.g.,
$Z_{5}~{}\Rightarrow{}~{}u_{15}=2u_{16}$ and $\neg
Z_{5}~{}\Rightarrow{}~{}u_{15}=u_{16}=0$, are converted into linear big-M
constraints, with $M>0$ fixed and large enough. For disjunctions with
nonlinear expressions we introduce auxiliary variables $w$ and augment the
(otherwise linear) objective with a quadratic penalty term: for instance,
$Z_{1}~{}\Rightarrow{}~{}\varphi_{1}(u)\leq 0$ becomes
$w_{1}\leq{}M(1-z_{1})$, with the nonlinear term
$\lambda[\varphi_{1}(u)-w_{1}]^{2}$ added to the objective function, for some
large and fixed penalty parameter $\lambda>0$.

#### Setup

Consider the problem of [15, Example 2], reformulated as discussed above with
parameters $\lambda=M=10^{3}$; in the form (P) there are 30 real-valued and 8
binary-valued decision variables, 11 linear equality constraints, 37 linear
inequality constraints, and simple bounds. We run Algorithm 3.1 on this
instance starting from 100 initial points obtained by sampling a standard
normal distribution; we then execute the variant of Algorithm 3.1 with NLP
refinement from the same initial points. To stop the execution, we limit the
number of successful iterations (that is, actual updates) to 1000.

#### Results

Algorithm 3.1 does not solve any of the 100 runs to the specified tolerance,
hitting the iteration limit in 96 cases and otherwise generating a negative
criticality measure. In contrast, the variant with NLP refinement returns
successfully in 69 cases; the MILP solver declares the subproblem infeasible
in 8 cases and generates a negative criticality measure in the remaining 23
cases. The best solution (in terms of objective value) is visualized in Figure
3, while computational results are summarized in Table 2.

Without NLP refinement, the algorithm was not able to reach an approximate
critical point within the tolerance. Behaving like a steepest descent method
when such steps are skipped, Algorithm 3.1 suffers from the slow tail
convergence that is typical of first-order methods. In contrast, despite the
computational cost of additional NLP solves, the algorithm can converge
quickly when using second-order information. It is particularly important,
however, that the integer-valued variables identify soon and correctly the
active set. The relatively low number of iterations reported in Table 2
indicate the potential benefits of NLP refinement steps within Algorithm 3.1
to speed up convergence.

$z_{1}$$z_{2}$$z_{3}$$z_{4}$$z_{5}$$z_{6}$$z_{7}$$z_{8}$$u_{1}$$u_{2}$$u_{3}$$u_{4}$$u_{5}$$u_{11}$$u_{12}$$u_{15}$$u_{6}$$u_{7}$$u_{8}$$u_{20}$$u_{22}$$u_{23}$$u_{24}$$u_{14}$$u_{13}$$u_{19}$$u_{21}$$u_{16}$$u_{17}$$u_{25}$$u_{9}$$u_{10}$$u_{18}$ Figure 3: Superstructure for the processing network example: units (boxes) and edges (lines) selected (solid) and discarded (dashed) by the optimal design. |  | min | q(25%) | median | q(75%) | max
---|---|---|---|---|---|---
without NLP | objective value | 76.92 | 80.31 | 90.62 | 103.60 | 116.39
refinement | iterations | 353 | 1882 | 1892.5 | 1907.5 | 1933
| successful updates | 186 | 1000# | 1000# | 1000# | 1000#
| runtime [s] | 1.27 | 6.15 | 6.19 | 6.24 | 6.47
with NLP | objective value | 59.85 | 72.09 | 83.75 | 95.65 | 107.56
refinement | iterations | 6 | 32.5 | 45.5 | 54 | 73
| successful updates | 2 | 12.5 | 18 | 22 | 29
| runtime [s] | 0.05 | 0.23 | 0.31 | 0.39 | 0.76
|  $\operatorname{proj}(x^{0})$ runtime [s] | 0.007 | 0.009 | 0.010 | 0.012 | 0.088
Table 2: Summary of computational results for the problem of a processing
network. Statistics (minimum, maximum, and quartiles) over different initial
points. Algorithm 3.1’s runtime includes NLP solves but not the sharp
projection of infeasible initial points, denoted $\operatorname{proj}(x^{0})$.
The symbol # stands for a limit reached.

##  5 Final remarks

The results in this paper could be extended by integrating acceleration
schemes and exploring the effects of inexact subproblem solves on the
convergence. Weaker assumptions on the problem structure and stronger
optimality concepts could also be examined, for instance replacing
$\Psi(x;\Delta)$ with $\Psi(x;\Delta)/\Delta$ as criticality measure. Future
research may focus on active-set warm-starting and lazification techniques for
computational efficiency, taking advantage of the similarity of successive
inner MILPs. The approach could also be adopted to tackle problems involving
additional nonlinear constraints, by integrating it within sequential
partially-constrained optimization schemes, such as penalty and barrier
methods.

## References

  * [1] Aleksandr Y. Aravkin, Robert Baraldi, and Dominique Orban.  A proximal quasi-Newton trust-region method for nonsmooth regularized optimization.  SIAM Journal on Optimization, 32(2):900–929, 2022.
  * [2] Francesco Biral, Enrico Bertolazzi, and Paolo Bosetti.  Notes on numerical methods for solving optimal control problems.  IEEJ Journal of Industry Applications, 5(2):154–166, 2016.
  * [3] Alexander Bodard, Pieter Pas, and Panagiotis Patrinos.  PANTR: A proximal algorithm with trust-region updates for nonconvex constrained optimization.  IEEE Control Systems Letters, 7:2389–2394, 2023.
  * [4] Richard H. Byrd, Nicholas I. M. Gould, Jorge Nocedal, and Richard A. Waltz.  On the convergence of successive linear-quadratic programming algorithms.  SIAM Journal on Optimization, 16(2):471–489, 2005.
  * [5] Alberto De Marchi.  Proximal gradient methods beyond monotony.  Journal of Nonsmooth Analysis and Optimization, 4, 2023.
  * [6] Alberto De Marchi and Andreas Themelis.  Proximal gradient algorithms under local Lipschitz gradient continuity.  Journal of Optimization Theory and Applications, 194(3):771–794, 2022.
  * [7] Oliver Exler and Klaus Schittkowski.  A trust region SQP algorithm for mixed-integer nonlinear programming.  Optimization Letters, 1(3):269–280, 2007.
  * [8] Michael P. Friedlander, Nick I. M. Gould, Sven Leyffer, and Todd S. Munson.  A filter active-set trust-region method.  Technical Report ANL/MCS-P1456-0907, Argonne National Laboratory, 2007\.
  * [9] Ignacio E. Grossmann and Francisco Trespalacios.  Systematic modeling of discrete-continuous optimization models through generalized disjunctive programming.  AIChE Journal, 59(9):3276–3295, 2013.
  * [10] Deborah Hendrych, Hannah Troppens, Mathieu Besançon, and Sebastian Pokutta.  Convex mixed-integer optimization with Frank-Wolfe methods, 2023.  arXiv:2208.11010v5.
  * [11] Christian Kirches, Jeffrey Larson, Sven Leyffer, and Paul Manns.  Sequential linearization method for bound-constrained mathematical programs with complementarity constraints.  SIAM Journal on Optimization, 32(1):75–99, 2022.
  * [12] Armin Nurkanović and Moritz Diehl.  Continuous optimization for control of hybrid systems with hysteresis via time-freezing.  IEEE Control Systems Letters, 6:3182–3187, 2022.
  * [13] Rien Quirynen and Stefano Di Cairano.  Sequential quadratic programming algorithm for real-time mixed-integer nonlinear MPC.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 993–999. IEEE Press, 2021.
  * [14] Andreas Themelis, Lorenzo Stella, and Panagiotis Patrinos.  Forward-backward envelope for the sum of two nonconvex functions: Further properties and nonmonotone linesearch algorithms.  SIAM Journal on Optimization, 28(3):2274–2303, 2018.
  * [15] Metin Türkay and Ignacio E. Grossmann.  Logic-based MINLP algorithms for the optimal synthesis of process networks.  Computers & Chemical Engineering, 20(8):959–978, 1996.
  * [16] Ya-xiang Yuan.  Conditions for convergence of trust region algorithms for nonsmooth optimization.  Mathematical Programming, 31(2):220–228, 1985.


# [Re] Double Sampling Randomized Smoothing

Aryan Gupta
Department of Electrical Engineering
Indian Institute of Technology, Roorkee
<EMAIL_ADDRESS>
&Sarthak Gupta
Department of Mathematics
Indian Institute of Technology, Roorkee
<EMAIL_ADDRESS>
&Abhay Kumar
Department of Mathematics
Indian Institute of Technology, Roorkee
<EMAIL_ADDRESS>
&Harsh Dugar
Department of Mathematics
Indian Institute of Technology, Roorkee
<EMAIL_ADDRESS>

## Reproducibility Summary

### Scope of Reproducibility

This effort aims to reproduce and validate the experiments and robustness of
the Double Sampling Randomized Smoothing[1] certification technique introduced
in [1]. We consistently verify the improvement in robustness radius for the
DSRS certification compared to the standard Neyman-Pearson
Certification[2][3][4][5] as reported in [1]. We also study the effectiveness
of the certification through adversarial trained models, ablation studies and
new experiments.

### Methodology

We start the reproduction effort by using the code open-sourced by the
authors. We reproduce Tables 7 and 8 of the original paper. Further, we
refactor their code for evaluating the method’s effectiveness under various
circumstances. We implement new ablation studies to test the method and design
experiments to verify the claims made by the authors. We release our code as
open-source.

### Results

We reproduced the results of DSRS certification on Cifar-10 trained Resnet-110
model within 1.7% accuracy of the reported values, and about half times we
reproduced even higher accuracy than the reported ones in the paper. The
authors’ claim that DSRS provides consistently tighter robustness
certification than existing baselines, seems to be consistently verified. The
ablation studies help us to further strengthen the claims made by the authors
and get further insights.

### What was easy

The authors open-sourced the code for the paper. This made it easy to verify
the results and allowed us to easily make edits for ablation studies. The
mechanism was also mathematically well described in the paper, helping us
conceptually understand the methods.

### What was difficult

The primary difficulty that we faced was the limitation of our computational
resources as compared to the requirements for conducting the paper’s
experiments and further ablation studies. For instance, some of the
experiments took over 60 hours to run.

### Communication with original authors

The authors assisted us with key insights in this reproducibility challenge,
from providing us the main code of the paper to their supportive suggestions
to our ablation studies. We would like to acknowledge Prof. Bo Li, Linyi Li,
and others for their invaluable time. Their feedback for this report is
available in the summary 5.4.

##  1 Introduction

Verifying the robustness of deep learning models is an active area of
research. There have been several efforts towards this direction [6][7]. One
significant work was Certified Adversarial Robustness via Randomized
Smoothing[3]. While it provided a method for effectively verifying the robust
radius of a deep network, it was not able to scale with the increasing input
dimensions. To improve this, Double Sampling Randomized Smoothing provides a
mathematical framework which builds upon [3]. The paper proposes to use an
additional distribution to calculate a tighter robust radius around the
samples. This distribution provides additional data for improving the robust
radius and can certify radius $\Omega(\sqrt[]{d})$ under the $l_{2}$ norm.

##  2 Scope of reproducibility

Under this reproducibility effort, we focus on reproducing the original
results of the paper as well as conducting further experiments to validate the
mechanism in varying settings. The claims from the original paper, which we
attempt to verify, are listed as follows.

  * •

DSRS certifies $\Omega(\sqrt[]{d})$ under the $l_{2}$ norm being the input
dimension, under mild assumptions. This implies that DSRS is able to break the
curse of dimensionality of randomized smoothing(“curse of dimensionality” or
“$l_{\infty}$ barrier”)[2][8][9][10].

  * •

DSRS provides consistently tighter robustness certification than existing
baselines, including standard Neyman-Pearson based certification across
different models on MNIST and CIFAR-10.

Further, we try to evaluate the proposed method by conducting the following
experiments:

  * •

Testing the DSRS method by varying the amount of samples.

  * •

Evaluating the effects of adversarial training on the certified radius of the
smoothed classifier.

  * •

Evaluating different distributions for smoothing than the ones used in the
original work.

  * •

Observing the variations in variance of the distributions used to smooth the
classifier.

##  3 Methodology

For this reproducibility study, we begin with the open-sourced code provided
by the authors and refactor it to conduct ablations wherever necessary. To
verify the claims listed in Section 2 we make use of the MNIST [11] and
CIFAR-10 [12]. Initially, we recreate the experiments exactly as listed in [1]
and then conduct several ablation studies by varying different hyperparameters
and configurations in the proposed method.

###  3.1 DSRS Algorithm Description

Figure 1:  A high-level overview of the DSRS method. [1]

Double sampling randomized smoothing (DSRS) is a novel robustness
certification framework for randomized smoothing. Randomized smoothing (RS)
[3] adds noise to the input and uses majority voting among the model
predictions for multiple noised inputs to get the final label prediction.
Since the shifting of the center of noised input distribution does not change
the ranking of label predictions, the RS models are certifiably robust to
small $l_{2}$-norm-bounded perturbations.

Given an input instance and a model, robustness certification approach
computes a robust radius, such that any perturbations within the radius does
not the change the final prediction. By nature, certification approach is
conservative — it provides a lower bound of robust radius while the maximum
robust radius for the given input instance could be much larger than the one
computed by the certification approach.

For RS, the most widely-used certification approach is the Neyman-Pearson-
based certification[2][3][4][5]. It leverages the probability of base model’s
predicting each class under the input noise to compute the certification.

In DSRS, authors propose to sample the base model’s prediction statistics
under two different distributions, and leverage the joint information to
compute the certification. Since more information is being utilized, the
certification approach is guaranteed to be at least as tight as the most
widely-used Neyman-Pearson-based approach.

###  3.2 Hyperparameters

The main hyperparameters in the certification are related to the two sampling
distributions i.e. P and Q. We studied the effect of changing the
hyperparameter k of the generalized gaussian distributions[13]. Also we
conducted experiments on variances $\sigma_{P}$ and $\sigma_{Q}$ of the main
distributions and studied the effects of their variation on the results.

###  3.3 Experimental setup, Datasets and code

For our studies, we use the MNIST and CIFAR-10 datasets containing 60,000 and
50,000 train examples respectively. We created bash scripts for each setting
of the experiment depending on the datasets and variance of the sampling
distributions of P(0.25, 0.50, and 1.00) and Q(0.20, 0.40 and 0.80), and the
distributions for sampling noise i.e. Standard Gaussian[3][2], and Generalized
Gaussian[13](k=380 for Mnist and k=1530 for Cifar10).Finally we trained each
model on three commonly used or state-of-the-art training methods: Gaussian
augmentation [3], Consistency [14], and SmoothMix [15]. We run the main
experiments of the paper 3 times and provide 95% confidence intervals for the
results in the Section 4. Also we created bash scripts for the respective
ablation studies like training adversarial models and hyperparameter search to
allow easy reproduction of experiments.

###  3.4 Computational requirements

We use a single NVIDIA Tesla V100 GPU with 32 GBs of RAM for this study. The
computational requirements for validating the results and conducting ablations
were high when compared to the compute power at our disposal. The
certification of robust radius for MNIST took over ten hours in all cases. For
CIFAR-10, the certification took over twenty-five hours, with some experiments
taking over sixty-five hours.

##  4 Results

In this section, we present the reproduction of certification results and
tables associated with the MNIST and CIFAR-10 datasets that the original paper
used to support their method. We could not use the larger ImageNet dataset due
to limited computational resources. We cover the results for MNIST and
CIFAR-10 in the same form as given in the original paper and conduct some
ablations to further investigate the DSRS framework.

###  4.1 Results reproducing original paper

Reproducing the original results of the paper helps us understand the main
claims of the paper, as mentioned in Section 2, and the main advantages of
DSRS certification over Standard Neyman-Pearson certification.

####  4.1.1 Verification of robustness radius having order
$\Omega(\sqrt[]{d})$

The DSRS framework seems to work well on increasing dimensions and provides us
with an increased ACR(Average Certified Radius)[16] as compared to the
standard Neyman-Pearson certification. Also, we can see the effect of varying
sigma values on the ACR values.

Figure 2:  The main result of DSRS method.[1]

####  4.1.2 Results verifying higher certified radius for DSRS on MNIST and
CIFAR10

The major results of the main paper are reproduced in this section. We have
reproduced these results three times and have calculated the confidence
intervals appropriately. Also, insights have been made into the results and
their variation from the ones provided in the main paper.

Table 1: MNIST: General Gaussian with 95% confidence interval

Variance |  Training Method |  Certification Approach | Certified Accuracy under Radius r
---|---|---|---
0.25 | 0.50 | 0.75 | 1.00 | 1.25 | 1.50 | 1.75 | 2.00 | 2.25 | 2.50 | 2.75 | 3.00
0.25 |  Gaussian Augmentation | Neyman-Pearson | 97.9% | 96.4% | 92.1% |  |  |  |  |  |  |  |  |
DSRS (OP) | 97.9% | 96.6% | 92.7% |  |  |  |  |  |  |  |  |
DSRS(Ours) | 97.8±0.3% | 96.7±0.1%* | 93.0±0.4%* |  |  |  |  |  |  |  |  |
Consistency | Neyman-Pearson | 98.4% | 97.5% | 94.4% |  |  |  |  |  |  |  |  |
DSRS (OP) | 98.4% | 97.5% | 95.4% |  |  |  |  |  |  |  |  |
DSRS(Ours) | 98.3±0.1% | 97.5±0.3% | 95.1±0.3% |  |  |  |  |  |  |  |  |
SmoothMix | Neyman-Pearson | 98.6% | 97.6% | 96.5% |  |  |  |  |  |  |  |  |
DSRS (OP) | 98.6% | 97.7% | 96.8% |  |  |  |  |  |  |  |  |
DSRS(Ours) | 98.6±0.2% | 97.2±0.3% | 95.5±0.2% |  |  |  |  |  |  |  |  |
0.50 |  Gaussian Augmentation | Neyman-Pearson | 97.8% | 96.9% | 94.6% | 88.4% | 78.7% | 52.6% |  |  |  |  |  |
DSRS (OP) | 97.8% | 97.0% | 95.0% | 89.8% | 83.4% | 59.1% |  |  |  |  |  |
DSRS(Ours) | 97.9±0.2%* | 97.0±0.1% | 95.3±0.5%* | 90.5±0.2%* | 82.6±1.0% | 59.8±1.8%* |  |  |  |  |  |
Consistency | Neyman-Pearson | 98.4% | 97.3% | 96.0% | 92.3% | 83.8% | 67.5% |  |  |  |  |  |
DSRS (OP) | 98.4% | 97.3% | 96.0% | 93.5% | 87.1% | 71.8% |  |  |  |  |  |
DSRS(Ours) | 98.3±0.3% | 97.3±0.2% | 96.1±0.3%* | 93.4±0.4% | 87.0±0.4% | 71.2±0.8% |  |  |  |  |  |
SmoothMix | Neyman-Pearson | 98.2% | 97.1% | 95.4% | 91.9% | 85.1% | 73.0% |  |  |  |  |  |
DSRS (OP) | 98.1% | 97.1% | 95.9% | 93.4% | 87.5% | 76.6% |  |  |  |  |  |
DSRS(Ours) | 97.5±0.1% | 96.4±0.1% | 94.5±0.1% | 91.9±0.2% | 86.3±0.2% | 75.4±0.6% |  |  |  |  |  |
1.00 |  Gaussian Augmentation | Neyman-Pearson | 95.2% | 91.9% | 87.7% | 80.6% | 71.2% | 57.6% | 41.0% | 25.5% | 13.6% | 6.2% | 2.1% | 0.9%
DSRS (OP) | 95.1% | 91.8% | 88.2% | 81.5% | 73.6% | 61.6% | 48.4% | 34.1% | 21.0% | 10.6% | 4.4% | 1.2%
DSRS(Ours) | 94.8±0.9% | 91.2±0.9% | 86.8±1.7% | 80.3±1.6% | 72.5±3.3% | 60.6±2.6% | 47.9±1.8% | 32.7±2.2% | 19.6±2.7% | 9.7±1.5% | 4.4±0.2% | 1.1±0.2%
Consistency | Neyman-Pearson | 93.9% | 90.9% | 86.4% | 80.8% | 73.0% | 61.1% | 49.1% | 35.6% | 21.7% | 10.4% | 4.1% | 1.9%
DSRS (OP) | 93.9% | 91.1% | 86.9% | 81.7% | 75.2% | 65.6% | 55.8% | 41.9% | 31.4% | 17.8% | 8.6% | 2.8%
DSRS(Ours) | 94.2±0.2%* | 91.2±0.1%* | 87.3±0.2%* | 81.7±0.3% | 74.7±0.2% | 65.6±0.7% | 54.7±1.0% | 41.9±0.8% | 30.2±1.1% | 17.1±1.5% | 8.6±0.1% | 2.6±0.1%
SmoothMix | Neyman-Pearson | 92.0% | 88.9% | 84.4% | 78.6% | 69.8% | 60.7% | 49.9% | 40.2% | 31.5% | 22.2% | 12.2% | 4.9%
DSRS (OP) | 92.2% | 89.0% | 84.8% | 79.7% | 72.0% | 63.9% | 54.4% | 46.2% | 37.6% | 29.2% | 18.5% | 7.2%
DSRS(Ours) | 91.9±0.1% | 88.7±0.2% | 85.0±0.1% | 79.8±0.3%* | 72.4±0.2%* | 63.9±0.1% | 54.9±0.2%* | 46.1±0.3% | 37.5±0.1% | 29.0±0.2% | 19.1±0.3%* | 7.2±0.2%

Table 2: CIFAR-10: General Gaussian with 95% confidence interval

Variance |  Training Method |  Certification Approach | Certified Accuracy under Radius r
---|---|---|---
0.25 | 0.50 | 0.75 | 1.00 | 1.25 | 1.50 | 1.75 | 2.00 | 2.25 | 2.50 | 2.75 | 3.00
0.25 |  Gaussian Augmentation | Neyman-Pearson | 56.1% | 35.7% | 13.4% |  |  |  |  |  |  |  |  |
DSRS (OP) | 57.4% | 39.4% | 17.3% |  |  |  |  |  |  |  |  |
DSRS(Ours) | 57.7±0.8%* | 37.8±0.7% | 17.5±0.4%* |  |  |  |  |  |  |  |  |
Consistency | Neyman-Pearson | 61.8% | 50.9% | 34.7% |  |  |  |  |  |  |  |  |
DSRS (OP) | 62.5% | 52.5% | 37.8% |  |  |  |  |  |  |  |  |
DSRS(Ours) | 62.8±0.7%* | 52.9±0.8%* | 39.2±0.9%* |  |  |  |  |  |  |  |  |
SmoothMix | Neyman-Pearson | 63.9% | 53.3% | 38.4% |  |  |  |  |  |  |  |  |
DSRS (OP) | 64.7% | 55.5% | 41.1% |  |  |  |  |  |  |  |  |
DSRS(Ours) | 60.5±0.6% | 53.2±0.6% | 42.5±0.5%* |  |  |  |  |  |  |  |  |
0.50 |  Gaussian Augmentation | Neyman-Pearson | 53.7% | 41.3% | 27.7% | 17.1% | 9.1% | 2.8% |  |  |  |  |  |
DSRS (OP) | 54.1% | 42.7% | 30.6% | 20.3% | 12.6% | 4.0% |  |  |  |  |  |
DSRS(Ours) | 52.8±0.8% | 40.6±1.4% | 29.6±1.5% | 20.2±0.9% | 12.0±0.3% | 4.7±0.5%* |  |  |  |  |  |
Consistency | Neyman-Pearson | 49.2% | 43.9% | 38.0% | 32.3% | 23.8% | 18.1% |  |  |  |  |  |
DSRS (OP) | 49.6% | 44.1% | 38.7% | 35.2% | 28.1% | 19.7% |  |  |  |  |  |
DSRS(Ours) | 50.3±0.8%* | 44.8±0.2%* | 39.8±0.6%* | 34.2±0.3% | 27.7±0.5% | 21.1±0.1%* |  |  |  |  |  |
SmoothMix | Neyman-Pearson | 53.2% | 47.6% | 40.2% | 34.2% | 26.7% | 19.6% |  |  |  |  |  |
DSRS (OP) | 53.3% | 48.5% | 42.1% | 35.9% | 29.4% | 21.7% |  |  |  |  |  |
DSRS(Ours) | 52.7±1.6% | 47.4±0.8% | 41.5±1.2% | 35.8±0.8% | 28.8±0.5% | 20.3±0.5% |  |  |  |  |  |
1.00 |  Gaussian Augmentation | Neyman-Pearson | 40.2% | 32.6% | 24.7% | 18.9% | 14.9% | 10.2% | 7.5% | 4.1% | 2.0% | 0.7% | 0.1% | 0.1%
DSRS (OP) | 40.3% | 33.1% | 25.9% | 20.6% | 16.1% | 12.5% | 8.4% | 6.4% | 3.5% | 1.8% | 0.7% | 0.1%
DSRS(Ours) | 40.1±1.0% | 32.4±0.5% | 26.1±0.2%* | 20.6±0.4% | 16.5±0.6%* | 12.4±0.4% | 8.9±0.5%* | 6.2±0.3% | 4.1±0.5%* | 2.1±0.3%* | 0.8±0.2%* | 0.2±0.1%*
Consistency | Neyman-Pearson | 37.2% | 32.6% | 29.6% | 25.9% | 22.5% | 19.0% | 16.4% | 13.8% | 11.2% | 9.0% | 7.1% | 5.1%
DSRS (OP) | 37.1% | 32.5% | 29.8% | 27.1% | 23.5% | 20.9% | 17.6% | 15.3% | 13.1% | 10.9% | 8.9% | 6.5%
DSRS(Ours) | 36.8±0.2% | 33.1±0.5%* | 29.5±0.9% | 26.2±0.7% | 23.5±0.7% | 20.2±0.1% | 17.5±0.5% | 15.7±0.2%* | 13.4±0.3%* | 11.1±0.3%* | 8.9±0.1% | 6.7±0.4%*
SmoothMix | Neyman-Pearson | 43.2% | 39.5% | 33.9% | 29.1% | 24.0% | 20.4% | 17.0% | 13.9% | 10.3% | 7.8% | 4.9% | 2.3%
DSRS (OP) | 43.2% | 39.7% | 34.9% | 30.0% | 25.8% | 22.1% | 18.7% | 16.1% | 13.2% | 10.2% | 7.1% | 3.9%
DSRS(Ours) | 43.3±1.8%* | 38.6±1.0% | 34.1±0.9% | 29.6±0.7% | 25.0±0.3% | 21.8±0.4% | 18.3±0.3% | 15.3±0.2% | 13.1±0.4% | 10.2±0.3% | 7.0±0.5% | 3.9±0.3%

###  4.2 Results beyond original paper

According to the original paper’s authors, the DSRS framework provides the
best results on the General-Gaussian distributions. In this section, we’ll
look at the effects of distribution shifts, specifically through Standard
Gaussian. Also, we are going to observe the effects of the different training
methods as suggested in the main paper, like PGD(Projected Gradient Descent)
Training[17] and Interpolated Training[18] on the Neyman-Pearson and DSRS
certified radius and compare these results. Finally, we are going to conduct
experiments on various hyperparameters introduced by the certification.

####  4.2.1 Results on Standard Gaussian

The following results are based on Standard-Gaussian distribution[3][2], and
we can see that DSRS performs equivalently to or a bit less than Neyman-
Pearson Certification. This is further explained in the ablations for the k
hyperparameter.

Table 3: CIFAR-10: Standard Gaussian

Variance |  Training Method |  Certification Approach | Certified Accuracy under Radius r
---|---|---|---
0.25 | 0.50 | 0.75 | 1.00 | 1.25 | 1.50 | 1.75 | 2.00 | 2.25 | 2.50 | 2.75 | 3.00
0.25 |  Gaussian Augmentation | Neyman-Pearson | 57.9% | 40.2% | 24.2% |  |  |  |  |  |  |  |  |
DSRS | 57.9% | 39.7% | 23.4% |  |  |  |  |  |  |  |  |
Consistency | Neyman-Pearson | 63.0% | 53.9% | 42.6% |  |  |  |  |  |  |  |  |
DSRS | 63.0% | 53.7% | 42.4% |  |  |  |  |  |  |  |  |
SmoothMix | Neyman-Pearson | 61.6% | 54.5% | 46.7% |  |  |  |  |  |  |  |  |
DSRS | 61.5% | 54.5% | 46.2% |  |  |  |  |  |  |  |  |
0.50 |  Gaussian Augmentation | Neyman-Pearson | 51.2% | 41.2% | 31.5% | 20.8% | 13.3% | 8.8% |  |  |  |  |  |
DSRS | 51.0% | 41.1% | 31.4% | 20.4% | 13.2% | 8.5% |  |  |  |  |  |
Consistency | Neyman-Pearson | 49.8% | 44.9% | 39.8% | 35.1% | 31.4% | 25.2% |  |  |  |  |  |
DSRS | 49.7% | 44.7% | 39.5% | 35.0% | 30.1% | 24.8% |  |  |  |  |  |
SmoothMix | Neyman-Pearson | 54.5% | 47.5% | 41.9% | 35.8% | 30.0% | 25.5% |  |  |  |  |  |
DSRS | 54.4% | 47.5% | 41.7% | 35.7% | 30.0% | 24.8% |  |  |  |  |  |
1.00 |  Gaussian Augmentation | Neyman-Pearson | 38.4% | 32.7% | 25.9% | 20.8% | 16.7% | 12.9% | 9.3% | 7.2% | 4.6% | 3.2% | 2.1% | 1.0%
DSRS | 38.3% | 32.4% | 25.8% | 20.7% | 16.7% | 12.7% | 9.0% | 7.0% | 4.5% | 3.2% | 2.0% | 0.9%
Consistency | Neyman-Pearson | 37.3% | 33.1% | 29.8% | 26.3% | 23.6% | 21.0% | 18.2% | 16.1% | 14.5% | 12.6% | 10.6% | 8.9%
DSRS | 37.1% | 32.9% | 29.8% | 6.3% | 23.6% | 20.7% | 17.9% | 16.0% | 14.0% | 12.4% | 10.2% | 8.6%
SmoothMix | Neyman-Pearson | 40.9% | 36.2% | 31.9% | 28.0% | 24.5% | 21.5% | 18.4% | 15.5% | 13.8% | 11.5% | 9.0% | 7.4%
DSRS | 40.9% | 36.2% | 31.8% | 27.8% | 24.4% | 21.3% | 18.3% | 15.5% | 13.8% | 11.3% | 9.0% | 7.0%

####  4.2.2 Results of adversarial models on General Gaussian and Standard
Gaussian

We then conducted experiments on two, adversarial trained models, i.e.,
PGD(Projected Gradient Descent) Training[17] and Interpolated Training[18]
trained models and then find the Neyman-Pearson and DSRS-certified radius for
both of them.

For the PGD model, we trained Resnet-110[19] on PGD based training and
obtained a benign test accuracy of 82.72% and adversarial test accuracy of
49.12% in the initial training. For Interpolated training, we trained the same
for Resnet-110[19] and obtained a benign test accuracy of 89.55% and
adversarial test accuracy of 47.86% in the initial training. Then we trained
the same on the 0.5 noise using all three training methods, and then
calculated the certification results for both the Standard and General
Gaussian distributions .

Table 4: CIFAR-10: Adversarially trained Models on General Gaussian

Training Method | Model |  Certification Approach | Certified Accuracy under Radius r |  Average Certified Radius
---|---|---|---|---
0.25 | 0.50 | 0.75 | 1.00 | 1.25 | 1.50
Gaussian Augmentation | Normal | Neyman-Pearson | 50.7% | 38.8% | 26.8% | 16.4% | 8.2% | 2.7% | 0.435
DSRS | 52.8% | 40.6% | 29.% | 20.2% | 12.0% | 4.7% | 0.467
Interpolated | Neyman-Pearson | 59.1% | 45.4% | 32.3% | 20.8% | 10.5% | 5.3% | 0.523
DSRS | 59.5% | 47.4% | 35.1% | 25.6% | 15.7% | 7.4% | 0.547
PGD | Neyman-Pearson | 58.4% | 46.4% | 31.8% | 20.3% | 11.0% | 4.2% | 0.514
DSRS | 59.0% | 47.6% | 35.2% | 24.0% | 14.1% | 6.4% | 0.538
Consistency | Normal | Neyman-Pearson | 49.4% | 44.4% | 39.5% | 31.6% | 24.8% | 18.9% | 0.608
DSRS | 50.3% | 44.8% | 39.8% | 34.2% | 27.7% | 21.1% | 0.610
Interpolated | Neyman-Pearson | 56.4% | 49.7% | 42.3% | 34.4% | 26.0% | 18.3% | 0.659
DSRS | 56.4% | 50.5% | 44.2% | 37.7% | 29.4% | 20.4% | 0.669
PGD | Neyman-Pearson | 56.0% | 49.4% | 42.1% | 33.8% | 26.2% | 17.6% | 0.652
DSRS | 56.1% | 50.3% | 43.0% | 36.9% | 29.8% | 20.6% | 0.663
SmoothMix | Normal | Neyman-Pearson | 53.4% | 47.3% | 40.5% | 33.0% | 25.3% | 17.7% | 0.638
DSRS | 52.7% | 47.4% | 41.5% | 35.8% | 28.8% | 20.3% | 0.638
Interpolated | Neyman-Pearson | 55.7% | 48.9% | 41.3% | 34.4% | 27.3% | 19.1% | 0.662
DSRS | 56.3% | 50.0% | 43.3% | 37.4% | 30.2% | 21.5% | 0.670
PGD | Neyman-Pearson | 55.4% | 48.7% | 40.6% | 33.4% | 26.7% | 18.7% | 0.655
DSRS | 55.6% | 49.6% | 42.8% | 35.8% | 29.8% | 21.9% | 0.661

Table 5: CIFAR-10: Adversarially trained Models on Standard Gaussian

Training Method | Model |  Certification Approach | Certified Accuracy under Radius r |  Average Certified Radius
---|---|---|---|---
0.25 | 0.50 | 0.75 | 1.00 | 1.25 | 1.50
Gaussian Augmentation | Normal | Neyman-Pearson | 51.2% | 41.2% | 31.5% | 20.8% | 13.3% | 8.8% | 0.503
DSRS | 51.0% | 41.1% | 31.4% | 20.4% | 13.2% | 8.5% | 0.496
Interpolated | Neyman-Pearson | 58.8% | 47.4% | 35.6% | 26.2% | 17.7% | 11.2% | 0.594
DSRS | 58.6% | 47.2% | 35.4% | 26.2% | 17.6% | 10.0% | 0.583
PGD | Neyman-Pearson | 57.9% | 47.4% | 35.3% | 25.8% | 16.9% | 10.8% | 0.582
DSRS | 57.7% | 47.1% | 35.3% | 25.6% | 16.4% | 10.3% | 0.571
Consistency | Normal | Neyman-Pearson | 49.8% | 44.9% | 39.8% | 35.1% | 30.4% | 25.2% | 0.681
DSRS | 49.7% | 44.7% | 39.5% | 35.0% | 30.1% | 24.8% | 0.658
Interpolated | Neyman-Pearson | 56.4% | 50.4% | 44.3% | 38.3% | 32.2% | 25.7% | 0.743
DSRS | 56.4% | 50.3% | 44.2% | 38.1% | 32.0% | 25.3% | 0.720
PGD | Neyman-Pearson | 56.0% | 49.9% | 43.3% | 37.5% | 31.4% | 25.7% | 0.735
DSRS | 55.9% | 49.9% | 43.0% | 37.5% | 31.3% | 25.1% | 0.714
SmoothMix | Normal | Neyman-Pearson | 54.5% | 47.5% | 41.9% | 35.8% | 30.0% | 25.5% | 0.713
DSRS | 54.4% | 47.5% | 41.7% | 35.7% | 30.0% | 24.8% | 0.691
Interpolated | Neyman-Pearson | 55.6% | 50.0% | 43.8% | 37.9% | 32.3% | 25.9% | 0.742
DSRS | 55.5% | 49.7% | 43.6% | 37.8% | 32.1% | 25.2% | 0.719
PGD | Neyman-Pearson | 55.0% | 50.0% | 43.1% | 36.2% | 31.5% | 26.6% | 0.733
DSRS | 55.0% | 49.9% | 43.0% | 36.1% | 31.3% | 25.9% | 0.710

We found that the results for General Gaussian were far better for DSRS
Certification than the Neyman-Pearson but the same can’t be seen on Standard
Gaussian results. This suggests that the adversarial training methods can
influence the robustness radius to a great extent and thus these results bring
an interesting future research scope in this area. The results of Standard
Gaussian are due to the theoretical limitations of the DSRS framework and can
be explained in the ablations of the k hyperparameter, as shown below.

####  4.2.3 Effect of ACR on varying hyperparameters N and k

Figure 3:  Ablations on the main hyperparameters- N and k shows the insights
of the DSRS famework

Here we can see that for N ablations, the results seem to increase linearly
for low values of N, i.e., from 50 to 10000, but after that the ACR nearly
becomes constant and does not show an appreciable increase. This is in line
with similar effects observed on Neyman-Pearson certified radius, but we could
have expected DSRS radius to increase more for larger values of N given its
theoretical advantages for two sampling distributions.

Also, we can see the ablation on the hyperparameter k which indicates the
$\sigma^{\prime}=\sqrt{\frac{d}{d-2k}}\sigma$ of General Gaussian
Distributions[13]. Here, we can see that as we move from $0$ to $d/2-15$ the
Neyman-Pearson Certification consistently outperforms the DSRS radius, but it
appears to work the other way for the range $d/2-15$ to $d/2$, as theorized by
the authors. As a result, the DSRS framework is better suited for higher
variance-based gaussian distributions as the parameter k changes the
$\sigma^{\prime}$ which is the variance for generalised gaussian
distributions[1].

####  4.2.4 Effect of ACR on varying Sigma P and Sigma Q

We did experiments to study the distributional shift on $\sigma_{P}$ and
$\sigma_{Q}$ values. We saw similar results as the authors did with
$\sigma_{Q}=0.4$, $(\sigma_{P}=0.5)$, ACR comes out to be maximum.

Also for changing $\sigma_{P}$ and varying $\sigma_{Q}$ appropriately, we
observed the maximum radius for $\sigma_{P}=1.0$ which shows a direct relation
between higher sigma values of the P and Q distributions and the higher ACR
values.

Figure 4:  Variation of ACR for $\sigma$ of P and Q distributions

##  5 Discussion

The results obtained by us largely support the authors’ claims that the DSRS
verification method is able to provide a tighter robustness certification.
While we could not compare and verify the superiority in all cases due to
computational limitations, we found the proposed method effective in the
various settings we tested. The experiments evaluating the certified radius
for various amounts of samplings provide ample evidence that the mechanism can
provide comparable results even with much lower sampling numbers, making it
much more feasible for real-life deployment.

###  5.1 What was easy

The authors open-sourced the code for the paper. This made it easy to verify
many results reported in the paper. The code was clear and concise, allowing
us to easily make edits for ablation studies. The mechanism was also
mathematically well described in the paper. This helped us conceptually
understand the methods.

###  5.2 What was difficult

The primary difficulty was the excessive requirement of computing capacity for
the experiments. Each verification on the MNIST dataset with 50,000 samples
took around 10 hours while each verification on the CIFAR-10 dataset took over
25 hours on an NVIDIA V100 GPU. This restricted our selection of datasets and
the scope of experiments we were able to conduct.

###  5.3 Recommendations for Reproducibility

The authors provided very clear instructions for reproducing the main results
in the GitHub repository of the paper. They also included the python scripts
that can be used for reproducing the main ablations of the paper. They were
used to produce the Main Robustness Verification curve in Section 4.1.1. Also,
the code was well factored which can make the reproduction effort seemingly
easy. The only thing that could have been done better is to provide clear
instructions about the ablation scripts as we had a hard time debugging the
files for our use.

###  5.4 Communication with original authors

We had extensive communication with the original authors throughout the course
of the main reproduction of results as well as the ablations. One of the
authors, Linyi Li responded well to the PGD[17] and Interpolated training[18]
methods we introduced in this paper. Also he helped us with other major
ablation experiments and suggested some more experiments, one of which we have
included in the final results.

To make the assessment fair and complete, we emailed to the main authors of
the paper, particularly Prof. Bo Li and Linyi Li, for their kind feedback on
the main report. They gave a very positive feedback for the report and
appreciated the ablations we did on the training methods and the
hyperparameters. They also acknowledged our efforts and said, "this report, as
a role model, is an awesome add-on to our work and greatly benefits the ML
community".

## References

  * [1] L. Li, J. Zhang, T. Xie, and B. Li, “Double sampling randomized smoothing,” 2022\.
  * [2] G. Yang, T. Duan, J. E. Hu, H. Salman, I. Razenshteyn, and J. Li, “Randomized smoothing of all shapes and sizes,” 2020.
  * [3] J. M. Cohen, E. Rosenfeld, and J. Z. Kolter, “Certified adversarial robustness via randomized smoothing,” 2019.
  * [4] H. Salman, G. Yang, J. Li, P. Zhang, H. Zhang, I. Razenshteyn, and S. Bubeck, “Provably robust deep learning via adversarially trained smoothed classifiers,” 2019.
  * [5] B. Li, C. Chen, W. Wang, and L. Carin, “Certified adversarial robustness with additive noise,” 2018.
  * [6] J. Z. Kolter and E. Wong, “Provable defenses against adversarial examples via the convex outer adversarial polytope,” CoRR, vol. abs/1711.00851, 2017\.
  * [7] M. Lecuyer, V. Atlidakis, R. Geambasu, D. Hsu, and S. Jana, “Certified robustness to adversarial examples with differential privacy,” 2018.
  * [8] A. Blum, T. Dick, N. Manoj, and H. Zhang, “Random smoothing might be unable to certify $\ell_{\infty}$ robustness for high-dimensional images,” 2020.
  * [9] A. Kumar, A. Levine, T. Goldstein, and S. Feizi, “Curse of dimensionality on randomized smoothing for certifiable robustness,” 2020.
  * [10] J. Hayes, “Extensions and limitations of randomized smoothing for robustness guarantees,” 2020.
  * [11] L. Deng, “The mnist database of handwritten digit images for machine learning research,” IEEE Signal Processing Magazine, vol. 29, no. 6, pp. 141–142, 2012.
  * [12] A. Krizhevsky, V. Nair, and G. Hinton, “Cifar-10 (canadian institute for advanced research),”
  * [13] D. Zhang, M. Ye, C. Gong, Z. Zhu, and Q. Liu, “Black-box certification with randomized smoothing: A functional optimization based framework,” 2020.
  * [14] J. Jeong and J. Shin, “Consistency regularization for certified robustness of smoothed classifiers,” 2020.
  * [15] J. Jeong, S. Park, M. Kim, H.-C. Lee, D. Kim, and J. Shin, “Smoothmix: Training confidence-calibrated smoothed classifiers for certified robustness,” 2021.
  * [16] R. Zhai, C. Dan, D. He, H. Zhang, B. Gong, P. Ravikumar, C.-J. Hsieh, and L. Wang, “Macer: Attack-free and scalable robust training via maximizing certified radius,” 2020.
  * [17] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu, “Towards deep learning models resistant to adversarial attacks,” 2017.
  * [18] A. Lamb, V. Verma, K. Kawaguchi, A. Matyasko, S. Khosla, J. Kannala, and Y. Bengio, “Interpolated adversarial training: Achieving robust neural networks without sacrificing too much accuracy,” 2019.
  * [19] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” 2015.