# A Primer on Multi-Neuron Relaxation-based Adversarial Robustness
Certification

Kevin Roth
Department of Computer Science
ETH Zürich
<EMAIL_ADDRESS>
###### Abstract

The existence of adversarial examples poses a real danger when deep neural
networks are deployed in the real world. The go-to strategy to quantify this
vulnerability is to evaluate the model against specific attack algorithms.
This approach is however inherently limited, as it says little about the
robustness of the model against more powerful attacks not included in the
evaluation. We develop a unified mathematical framework to describe
relaxation-based robustness certification methods, which go beyond adversary-
specific robustness evaluation and instead provide provable robustness
guarantees against attacks by any adversary. We discuss the fundamental
limitations posed by single-neuron relaxations and show how the recent
“k-ReLU” multi-neuron relaxation framework of Singh et al., 2019a  obtains
tighter correlation-aware activation bounds by leveraging additional
relational constraints among groups of neurons. Specifically, we show how
additional pre-activation bounds can be mapped to corresponding post-
activation bounds and how they can in turn be used to obtain tighter
robustness certificates. We also present an intuitive way to visualize
different relaxation-based certification methods. By approximating multiple
non-linearities jointly instead of separately, the k-ReLU method is able to
bypass the convex barrier imposed by single neuron relaxations.

##  1 Introduction

#### Adversarial Examples

While deep neural networks have been used with great success for perceptual
tasks such as image classification or speech recognition, their performance
can deteriorate dramatically in the face of so-called adversarial examples
(Biggio et al.,, 2013; Szegedy et al.,, 2013; Goodfellow et al.,, 2014), i.e.
small specifically crafted perturbations of the input signal, often
imperceptible to humans, that are sufficient to induce large changes in the
model output, cf. Figure 1. Ever since their discovery, there has been a huge
interest in the machine learning community to try to understand their origin
and to mitigate their consequences.

|  |
---|---|---
Figure 1: Adversarial Examples: (left) original image “giant panda” (ImageNet
class 339), (middle) adversarial perturbation, (right) adversarially perturbed
image causing misclassification “goldfish” (ImageNet class 2) on a MobileNetV2
pretrained on ImageNet.

The most direct and popular strategy of robustification, called adversarial
training, aims to improve the robustness of a model by training it against an
adversary that perturbs each training example before passing it to the model
(Goodfellow et al.,, 2014; Kurakin et al.,, 2016; Madry et al.,, 2018). Only
recently it was shown that projected gradient ascent based adversarial
training is equivalent to data-dependent operator norm regularization (Roth et
al.,, 2020), confirming the long-standing argument that the sensitivity of a
neural network to adversarial examples is tied to its spectral properties.

A different strategy of defense is to detect whether or not the input has been
perturbed, by detecting characteristic regularities either in the adversarial
perturbations themselves or in the network activations they induce (Grosse et
al.,, 2017; Feinman et al.,, 2017; Xu et al.,, 2017; Metzen et al.,, 2017;
Carlini and Wagner, 2017a, ; Roth et al.,, 2019).

#### Robustness Certification

The existence of adversarial examples poses a real danger when deep neural
networks are deployed in the real world. The go-to strategy to quantify this
vulnerability is to evaluate the model against specific attack algorithms.
This approach is however inherently limited, as it says little about the
robustness of the model against more powerful attacks not included in the
evaluation (Carlini and Wagner, 2017b, ; Athalye et al.,, 2018). We therefore
need to go beyond adversary-specific robustness evaluation and instead provide
provable robustness guarantees against attacks by any adversary.

The idea of robustness certification is to find the largest neighborhood
(typically norm-bounded) that guarantees that no perturbation inside it can
change the network’s prediction, or equivalently to find the minimum
distortion required to induce a misprediction. Unfortunately, exactly
certifying the robustness of a network is an NP-complete problem (Katz et
al.,, 2017; Weng et al.,, 2018). Consider a ReLU network: to find the exact
mimimum distortion, two branches have to be considered for each ReLU
activation where the input can take both positive and negative values. This
makes exact verification methods computationally demanding even for small
networks (Katz et al.,, 2017; Weng et al.,, 2018; Zhang et al.,, 2018).

#### Exact certification

Exact combinatorial verifiers solve the robustness certification problem for
(piecewise linear) ReLU networks by employing either mixed integer linear
programming (MILP) solvers (Cheng et al.,, 2017; Lomuscio and Maganti,, 2017;
Tjeng et al.,, 2018), which use binary variables to encode the states of ReLU
activations, or satisfiability modulo theories (SMT) solvers (Katz et al.,,
2017; Carlini et al.,, 2017; Ehlers,, 2017), which encode the network into a
set of linear constraint satisfaction problems based on the ReLU branches.
Exact approaches guarantee to find the minimum distortion, however, due to the
NP-completeness of the problem, they are fundamentally limited to small
networks. For instance, it can take Reluplex (Katz et al.,, 2017) several
hours to verify a feedforward network with 5 inputs, 5 outputs and 300 hidden
neurons on a single data point.

#### Relaxation-based certification

While exact certification is hard, providing a guaranteed certified lower
bound for the minimum adversarial perturbation, resp. an upper bound on the
robust error (i.e. the probability of misprediction for a given uncertainty
set), can be done efficiently with relaxation-based certification methods
(Hein and Andriushchenko,, 2017; Wong and Kolter,, 2018; Raghunathan et al.,,
2018; Dvijotham et al.,, 2018; Weng et al.,, 2018; Zhang et al.,, 2018; Singh
et al.,, 2018; Mirman et al.,, 2018; Singh et al., 2019b, ; Qin et al.,, 2019;
Salman et al.,, 2019). Relaxation-based verifiers trade off precision
(increased false negative rates) with efficiency and scalability by convexly
relaxing (over-approximating) the non-linearities in the network. The
resulting certificates are still sound: they may fail to verify robustness for
a data point that is actually robust, but they never falsely certify a data
point that is not robust.

#### Single-Neuron Relaxation Barrier

The effectiveness of existing single-neuron relaxation-based verifiers is
inherently limited by the optimal convex relaxation obtainable by processing
each non-linearity separately (Salman et al.,, 2019). Relaxing the non-
linearities separately comes at the cost of losing correlations between
neurons: While existing frameworks do consider correlations between units to
compute bounds for higher layers, the bounds for neurons in a given layer are
computed individually, i.e. without interactions within the same layer (as a
result of considering separate worst-case perturbations for each neuron)
(Salman et al.,, 2019). Moreover, as the activation bounds are obtained
recursively, there is a risk that the error amplifies across layers, which is
particularly problematic for deep neural networks. For ReLU non-linearities,
we incur an over-approximation error each time an unstable neuron (taking both
positive and negative values) is relaxed.

#### Multi-Neuron Relaxation

The single-neuron relaxation barrier can be bypassed by considering multi-
neuron relaxations, i.e. by approximating multiple non-linearities jointly
instead of separately, as suggested in the recent LP-solver based “k-ReLU”
relaxation framework of Singh et al., 2019a . The k-ReLU relaxation framework
obtains tighter activation bounds by leveraging additional relational
constraints among groups of neurons, enabling significantly more precise
certification than existing state-of-the-art verifiers while still maintaining
scalability (Singh et al., 2019a, ). Note that the k-ReLU method has recently
been generalized to work for arbitrary activation functions (Müller et al.,,
2021).

Finally, we would like to mention that there is another approach to bypass the
single-neuron relaxation barrier: (Tjandraatmadja et al.,, 2020) obtain
tightened single neuron relaxations by considering bounds in terms of the
multivariate post-activations preceding the affine layer instead of the
univariate input to the ReLU non-linearity. However, since their approach is
conceptually quite different from the multi-neuron relaxation approach
discussed here, we will not discuss it further.

#### Robust training

Robust training procedures aim to integrate verification into the training of
neural networks so they are guaranteed to exhibit a certain level of
robustness to norm-bounded perturbations on the training set. As long as the
certification objective is differentiable with respect to the network
parameters, it can in principle be used to train verifiable neural networks
(Wong and Kolter,, 2018; Wong et al.,, 2018). Differentiating through
certificates is one way of robustifying the network during training, another
is adversarial training or data-dependent operator norm regularization. Some
of the key questions are to what extent the robustness certification
generalizes from verified to unverified examples and to what extend
adversarial training resp. data-dependent operator norm regularization improve
the robustness certifiability of deep neural networks.

##  2 Single-Neuron Relaxation

#### Notation.

Let $f:\mathcal{X}\subset{\mathbb{R}}^{n_{0}}\to{\mathbb{R}}^{n_{L}}$ be an
$L$-layer deep feedforward neural network, given by the equations111This
formulation captures all network architectures in which units in one layer
receive inputs from units in the previous layer, including fully-connected and
convolutional networks. To capture residual networks, in which units receive
inputs from units in multiple previous layers, the right hand side of the
equation would have to be extended to include dependencies on
${\mathbf{x}}^{(j)}$ for $j<i-1$.

$\displaystyle{\mathbf{x}}^{(i)}={\mathbf{W}}^{(i)}\sigma({\mathbf{x}}^{(i-1)})+{\mathbf{b}}^{(i)}\quad\textnormal{for}\quad
i\in\llbracket L\rrbracket\,,$ (1)

with input ${\mathbf{x}}^{(0)}={\mathbf{x}}$ and output
${\mathbf{x}}^{(L)}\equiv f({\mathbf{x}})$, where
${\mathbf{W}}^{(i)},{\mathbf{b}}^{(i)}$ denote the layer-wise weight matrix
and bias vector, and where $\sigma(\cdot)$ denotes the non-linear activation
function, with the convention that $\sigma({\mathbf{x}}^{(0)})={\mathbf{x}}$
is the identity activation. We use
$n_{i}\coloneqq\textnormal{dim}({\mathbf{x}}^{(i)})$ to denote the number of
neurons in layer $i$, and ${\mathbf{x}}^{\llbracket
L_{0},L_{1}\rrbracket}\coloneqq\\{{\mathbf{x}}^{(L_{0})},\dots,{\mathbf{x}}^{(L_{1})}\\}$
to denote the collection of all ${\mathbf{x}}^{(i)}$ for $i\in\llbracket
L_{0},L_{1}\rrbracket$, where $\llbracket L_{0},L_{1}\rrbracket$ denotes the
set of indices $\llbracket
L_{0},L_{1}\rrbracket\coloneqq\\{L_{0},\dots,L_{1}\\}$, with the shorthand
notation $\llbracket L\rrbracket\equiv\llbracket 1,L\rrbracket$. Note that in
our notation the pre-activations ${\mathbf{x}}^{(i)}$ have the same index as
the weight matrix ${\mathbf{W}}^{(i)}$ and bias vector ${\mathbf{b}}^{(i)}$ on
which they directly depend. Together, the above equations define the mapping

$\displaystyle f(\cdot)=$
$\displaystyle{\mathbf{W}}^{(L)}\sigma(\,\cdots\,{\mathbf{W}}^{(2)}\sigma({\mathbf{W}}^{(1)}(\cdot)+{\mathbf{b}}^{(1)})+{\mathbf{b}}^{(2)}\,\cdots\,)+{\mathbf{b}}^{(L)}\,.$
(2)

We also introduce the following useful notation:

######  Definition 1.

(Row / column indexing). Let $\mathbf{A}\in{\mathbb{R}}^{m\times n}$ be an
arbitrary real-valued matrix. We use $\mathbf{A}_{i\cdot\cdot}$ to denote the
$i$-th row and $\mathbf{A}_{:j}$ to denote the $j$-th column of $\mathbf{A}$,
for $i\in\llbracket m\rrbracket$ and $j\in\llbracket n\rrbracket$.

######  Definition 2.

(Positive / negative entries). Let $\mathbf{A}\in{\mathbb{R}}^{m\times n}$ be
an arbitrary real-valued matrix. Define

$\displaystyle\left[\,\cdot\,\right]_{+}:\mathbf{A}\to\mathbf{A}_{+}\coloneqq\max(\mathbf{A},0)\,,$
(3)
$\displaystyle\left[\,\cdot\,\right]_{-}:\mathbf{A}\to\mathbf{A}_{-}\coloneqq\min(\mathbf{A},0)\,,$

where the $\max(\cdot,\cdot)$ and $\min(\cdot,\cdot)$ are taken _entrywise_.
Note that by definition, $\mathbf{A}=\mathbf{A}_{+}+\mathbf{A}_{-}$.

#### Robustness Certification.

The network is considered certifiably robust with respect to input
${\mathbf{x}}$ and uncertainty set $\mathcal{B}({\mathbf{x}})$ if there is no
perturbation within $\mathcal{B}({\mathbf{x}})$ that can change the network’s
prediction, i.e. if the largest logit $f_{\hat{k}({\mathbf{x}})}$ remains
larger than any other logit $f_{k}$ within the entire uncertainty set
$\mathcal{B}({\mathbf{x}})$. Formally, the network is certifiably robust at
${\mathbf{x}}$ with respect to $\mathcal{B}({\mathbf{x}})$ if

$\displaystyle\min_{{\mathbf{x}}^{*}\in\mathcal{B}({\mathbf{x}})}f_{\hat{k}({\mathbf{x}})}({\mathbf{x}}^{*})-f_{k}({\mathbf{x}}^{*})>0\,,\quad\forall\,k\neq\hat{k}({\mathbf{x}})\,,$
(4)

where $\hat{k}({\mathbf{x}})=\operatorname*{arg\,max}_{j}f_{j}({\mathbf{x}})$.

Typical choices for $\mathcal{B}({\mathbf{x}})$ are the $\ell_{p}$-norm balls
$\mathcal{B}_{p}({\mathbf{x}};\epsilon)\coloneqq\\{{\mathbf{x}}^{*}:||{\mathbf{x}}^{*}-{\mathbf{x}}||_{p}\leq\epsilon\\}$
or the non-uniform box domain
$\mathcal{B}({\mathbf{x}};\\{\epsilon^{-}_{j},\epsilon^{+}_{j}\\}_{j\in\llbracket
n_{0}\rrbracket})\coloneqq\\{{\mathbf{x}}^{*}:{\mathbf{x}}_{j}-\epsilon^{-}_{j}\leq{\mathbf{x}}^{*}_{j}\leq{\mathbf{x}}_{j}+\epsilon^{+}_{j}\\}$
with $\epsilon^{-}_{j},\epsilon^{+}_{j}\in{\mathbb{R}}_{0}^{+}$.

The optimization domain in the above certification problem can be narrowed
down significantly if we are given lower- and upper- pre-activation bounds,
${\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}$,
for $i\in\llbracket L\\!-\\!1\rrbracket$. We will see below how (approximate)
pre-activation bounds can be computed. The corresponding optimization problem,
subsequently referred to via the short-hand notation
$\mathcal{O}(\mathbf{c},c_{0},L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$,
reads

$\displaystyle\min_{{\mathbf{x}}^{\llbracket 0,L\rrbracket}}$
$\displaystyle\mathbf{c}^{\top}{\mathbf{x}}^{(L)}+c_{0}$ (5) s.t.
$\displaystyle{\mathbf{x}}^{(i)}={\mathbf{W}}^{(i)}\sigma({\mathbf{x}}^{(i-1)})+{\mathbf{b}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket L\rrbracket$
$\displaystyle{\mathbf{x}}^{(0)}\in\mathcal{B}({\mathbf{x}})\,,\
{\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket L\\!-\\!1\rrbracket$

where the vector-inequalities are considered to hold element-wise, i.e.
${\boldsymbol{\ell}}^{(i)}_{j}\leq{\mathbf{x}}^{(i)}_{j}\leq{\boldsymbol{u}}^{(i)}_{j}$,
$\forall\,j\in\llbracket n_{i}\rrbracket$. For
$\mathbf{c}\in\\{\mathbf{e}_{\hat{k}({\mathbf{x}})}-\mathbf{e}_{k}\\}_{k\neq\hat{k}({\mathbf{x}})}$,
$c_{0}=0$ and
${\boldsymbol{\ell}}^{(i)}_{j}=-\infty,{\boldsymbol{u}}^{(i)}_{j}=\infty$,
this formulation is equivalent to Equation 4.

We denote the optimal value of
$\mathcal{O}(\mathbf{c},c_{0},L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$ by
$p^{*}_{\mathcal{O}}$. If $p^{*}_{\mathcal{O}}>0$ for all
$\mathbf{c}\in\\{\mathbf{e}_{\hat{k}({\mathbf{x}})}-\mathbf{e}_{k}\\}_{k\neq\hat{k}({\mathbf{x}})}$,
with $c_{0}=0$ and valid pre-activation bounds
${\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket}$, the
network is certifiably robust with respect to ${\mathbf{x}}$ and
$\mathcal{B}({\mathbf{x}})$.

Exact pre-activation bounds ${\boldsymbol{\ell}}^{(i)}_{j}$ and
${\boldsymbol{u}}^{(i)}_{j}$ are given by the minimization ($+\mathbf{e}_{j}$)
resp. maximization ($-\mathbf{e}_{j}$) problems
$\mathcal{O}(\pm\mathbf{e}_{j},0,i,{\boldsymbol{\ell}}^{\llbracket
i\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket i\\!-\\!1\rrbracket})$ for
all neurons $j\in\llbracket n_{i}\rrbracket$ and layers $i\in\llbracket
L\\!-\\!1\rrbracket$. Unfortunately, computing exact bounds is as NP-hard as
solving the exact certification problem.

Most of the existing robustness certification methods in use today are based
on a variant of one of two (pre)-activation bound computation paradigms: (i)
Interval Bound Propagation or (ii) Relaxation-based Bound Computation. The
pre-activation bounds play a crucial role in these certification methods: the
tighter they are, the lower the false-negative rate of the method.

#### Relaxation-based Certification

Exactly certifying the robustness of a network via Equation 5 is an NP-
complete problem (Katz et al.,, 2017; Weng et al.,, 2018), due to the non-
convex constraints imposed by the non-linearities. Relaxation-based verifiers
trade off precision with efficiency and scalability by convexly relaxing
(over-approximating) the non-linearities in the network.

The convex relaxation-based certification problem, subsequently referred to
via the short-hand notation
$\mathcal{C}(\mathbf{c},c_{0},L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$,
reads

$\displaystyle\min_{{\mathbf{x}}^{\llbracket
0,L\rrbracket},\,\,{\mathbf{z}}^{\llbracket 0,L\\!-\\!1\rrbracket}}$
$\displaystyle\mathbf{c}^{\top}{\mathbf{x}}^{(L)}+c_{0}$ (6) s.t.
$\displaystyle{\mathbf{x}}^{(i)}={\mathbf{W}}^{(i)}{\mathbf{z}}^{(i-1)}+{\mathbf{b}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket L\rrbracket$
$\displaystyle\underline{\sigma}({\mathbf{x}}^{(i)})\leq{\mathbf{z}}^{(i)}\leq\overline{\sigma}({\mathbf{x}}^{(i)})\,\
\textnormal{for}\ i\in\llbracket L\\!-\\!1\rrbracket$
$\displaystyle{\mathbf{x}}^{(0)}\in\mathcal{B}({\mathbf{x}})\,,\
{\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket L\\!-\\!1\rrbracket$

where $\underline{\sigma}({\mathbf{x}}^{(i)})$ resp.
$\overline{\sigma}({\mathbf{x}}^{(i)})$ are convex resp. concave bounding
functions satisfying
$\underline{\sigma}({\mathbf{x}}^{(i)})\leq\sigma({\mathbf{x}}^{(i)})\leq\overline{\sigma}({\mathbf{x}}^{(i)})$
for all
${\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}$,
and where ${\mathbf{z}}^{(i)}$ are the post-activation variables with the
convention that ${\mathbf{z}}^{(0)}={\mathbf{x}}^{(0)}$. We denote the optimal
value of $\mathcal{C}(\mathbf{c},c_{0},L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$ by
$p^{*}_{\mathcal{C}}$. Naturally, we have that $p^{*}_{\mathcal{C}}\leq
p^{*}_{\mathcal{O}}$. Thus, if $p^{*}_{\mathcal{C}}>0$ for all
$\mathbf{c}\in\\{\mathbf{e}_{\hat{k}({\mathbf{x}})}-\mathbf{e}_{k}\\}_{k\neq\hat{k}({\mathbf{x}})}$,
with $c_{0}=0$ and valid pre-activation bounds
${\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket}$, the
network is certifiably robust w.r.t. ${\mathbf{x}}$ and
$\mathcal{B}({\mathbf{x}})$.

#### Optimal Single-Neuron Relaxation

Under mild assumptions (non-interactivity, see below), the optimal convex
relaxation of a single non-linearity, i.e. its convex hull, is given by
$\underline{\sigma}_{\textnormal{opt}}(x)\leq\sigma(x)\leq\overline{\sigma}_{\textnormal{opt}}(x)$
(Salman et al.,, 2019), where

$\displaystyle\underline{\sigma}_{\textnormal{opt}}(x)\,\,\ \textnormal{is the
greatest convex function majored by $\sigma$}$ (7)
$\displaystyle\overline{\sigma}_{\textnormal{opt}}(x)\,\,\ \textnormal{is the
smallest concave function majoring $\sigma$}$

In general, for vector-valued non-linearities
$\sigma:\mathbb{R}^{n}\to\mathbb{R}^{m}$, the optimal convex relaxation may
not have a simple analytic form. However, if there is no interaction among the
outputs $\sigma_{i}$ for $i\in\llbracket m\rrbracket$, the optimal convex
relaxation does admit a simple analytic form (Salman et al.,, 2019):

######  Definition 3.

Salman et al.’s Definition B.2 (non-interactivity) Let
$\sigma:\mathbb{R}^{n}\to\mathbb{R}^{m}$ be a vector-valued non-linearity with
input
${\mathbf{x}}\in[{\boldsymbol{\ell}},{\boldsymbol{u}}]\subset\mathbb{R}^{n}$
and output $\sigma({\mathbf{x}})\in\mathbb{R}^{m}$. For each output
$\sigma_{j}({\mathbf{x}})$, let $I_{j}\subset\llbracket n\rrbracket$ be the
set of ${\mathbf{x}}$’s entries that affect $\sigma_{j}({\mathbf{x}})$. We
call the vector-valued non-linearity _non-interactive_ if the sets $I_{j}$ for
$j\in\llbracket m\rrbracket$ are mutually disjoint $I_{j}\cap I_{k}=\emptyset$
for all $j\neq k\in\llbracket m\rrbracket$.

All element-wise non-linearities such as (leaky)-ReLU, sigmoid and tanh are
non-interactive. MaxPooling is also non-interactive if the stride is no
smaller than the kernel size, i.e. if the receptive regions are non-
overlapping.

For the ReLU non-linearity $\sigma(x)=\max(x,0)$, the optimal convex
relaxation with respect to pre-activation bounds $\ell,u$, is given by the
triangle relaxation (Ehlers,, 2017),

$\displaystyle\underline{\sigma}(x)=\max(0,x)\
,\quad\overline{\sigma}(x)=\frac{u}{u-\ell}(x-\ell)$ (8)

see Figure 2 (middle) for an illustration.

We denote the optimal relaxation-based certification problem as
$\mathcal{C}_{\textnormal{opt}}$ and the corresponding optimal value of the
objective as $p^{*}_{\mathcal{C}_{\textnormal{opt}}}$.

#### Optimal LP-relaxed Verification

For piece-wise linear networks, including (leaky)-ReLU networks, the optimal
single-neuron relaxation-based certification problem
$\mathcal{C}_{\textrm{opt}}(\mathbf{c},c_{0},L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$ is a
linear programming problem and can thus be solved exactly with off-the-shelf
LP solvers. Two steps are required (Salman et al.,, 2019): (a) We first need
to obtain optimal single-neuron relaxation-based pre-activation bounds for all
neurons in the network except those in the logit layer. (b) We then solve the
LP-relaxed (primal or dual) certification problem exactly for the logit layer
of the network.

(a) Obtaining optimal single-neuron relaxation-based pre-activation bounds.
The optimal single-neuron relaxation-based pre-activation bounds are obtained
by recursively solving the minimization ($+\mathbf{e}_{j}$) resp. maximization
($-\mathbf{e}_{j}$) problems
$\mathcal{C}_{\textrm{opt}}(\pm\mathbf{e}_{j},0,i,{\boldsymbol{\ell}}^{\llbracket
i\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket i\\!-\\!1\rrbracket})$ for
all neurons $j\in\llbracket n_{i}\rrbracket$ at increasingly higher layers
$i\in\llbracket L\\!-\\!1\rrbracket$.

(b) Solving the LP-relaxed (primal or dual) certification problem for the
logit layer. We then solve the linear program
$\mathcal{C}_{\textrm{opt}}(\mathbf{c},0,L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$ for
all
$\mathbf{c}\in\\{\mathbf{e}_{\hat{k}({\mathbf{x}})}-\mathbf{e}_{k}\\}_{k\neq\hat{k}({\mathbf{x}})}$
with the above pre-activation bounds. If the solutions of all linear programs
are postive, i.e. if $p^{*}_{\mathcal{C}_{\textrm{opt}}}>0$ (primal) or
$d^{*}_{\mathcal{C}_{\textrm{opt}}}>0$ (dual) for all
$\mathbf{c}\in\\{\mathbf{e}_{\hat{k}({\mathbf{x}})}-\mathbf{e}_{k}\\}_{k\neq\hat{k}({\mathbf{x}})}$,
the network is certifiably robust w.r.t. ${\mathbf{x}}$ and
$\mathcal{B}({\mathbf{x}})$.

Note that the number of optimization sub-problems that need to be solved
scales linearly with the number of neurons, which can easily be in the
millions for deep neural networks (Salman et al.,, 2019). For this reason,
much of the literature on certification for deep neural networks has focused
on efficiently computing _approximate_ pre-activation bounds. An efficient and
scalable alternative to solving the optimal LP-relaxed verification problem
exactly is to only solve (b) exactly and instead of (a) to greedily compute
approximate but sound pre-activation bounds.

As we will see shortly, tighter pre-activation bounds also yield tighter
relaxations when over-approximating the non-linearities. Most of the existing
relaxation-based certification methods in use today are based on a variant of
one of two greedy (pre)-activation bound computation paradigms: (i) interval
bound propagation (Algorithm 1) or (ii) relaxation-based backsubstitution
(Algorithm 2). The pre-activation bounds play a crucial role in these
certification methods: the tighter they are, the lower the false-negative rate
of the method.

#### Interval Bound Propagation.

The simplest possible method to obtain approximate but sound pre-activation
bounds is given by the Interval Bound Propagation (IBP) algorithm (Dvijotham
et al.,, 2018), which is based on the idea that valid layer-wise pre-
activation bounds can be obtained by considering a separate worst-case
previous-layer perturbation ${\boldsymbol{\xi}}$ for each row
${\mathbf{W}}^{(i)}_{j\cdot\cdot}$, satisfying the constraint that
${\boldsymbol{\xi}}$ is within the previous-layer lower- and upper- pre-
activation bounds
${\boldsymbol{\ell}}^{(i-1)}_{j^{\prime}}\leq{\boldsymbol{\xi}}_{j^{\prime}}\leq{\boldsymbol{u}}^{(i-1)}_{j^{\prime}},\forall
j^{\prime}\in\llbracket n_{i-1}\rrbracket$. The $j^{\prime}$-th perturbation
entry of the greedy solution is uniquely determined by the sign of the
$j^{\prime}$-th entry of the vector ${\mathbf{W}}^{(i)}_{j\cdot\cdot}$. The
corresponding bounds are

$\displaystyle{\boldsymbol{\ell}}^{(i)}_{j}$
$\displaystyle\geq\big{[}{\mathbf{W}}^{(i)}_{j\cdot\cdot}\big{]}_{+}\sigma({\boldsymbol{\ell}}^{(i-1)})+\big{[}{\mathbf{W}}^{(i)}_{j\cdot\cdot}\big{]}_{-}\sigma({\boldsymbol{u}}^{(i-1)})+{\mathbf{b}}^{(i)}_{j}$
(9) $\displaystyle{\boldsymbol{u}}^{(i)}_{j}$
$\displaystyle\leq\big{[}{\mathbf{W}}^{(i)}_{j\cdot\cdot}\big{]}_{-}\sigma({\boldsymbol{\ell}}^{(i-1)})+\big{[}{\mathbf{W}}^{(i)}_{j\cdot\cdot}\big{]}_{+}\sigma({\boldsymbol{u}}^{(i-1)})+{\mathbf{b}}^{(i)}_{j}$

The computational complexity of the Interval Bound Propagation algorithm is
linear in the number of layers $L$. The complete algorithm is shown in Section
5.2 in the Appendix.

A more sophisticated class of approximate but sound pre-activation bounds can
be obtained by greedy relaxation-based backsubstitution.

#### Linear bounding functions

Of particular interest is the case where each non-linear layer is bounded by
exactly one linear lower relaxation function
$\underline{{\boldsymbol{\Phi}}}^{(i)}(\cdot)$ and one linear upper relaxation
function $\overline{{\boldsymbol{\Phi}}}^{(i)}(\cdot)$, as the corresponding
relaxation-based verification problem can be solved greedily in this case
(Salman et al.,, 2019). For ReLU non-linearities we have the following
expressions for the linear bounding functions:

######  Proposition 1.

(Linear ReLU relaxation functions (Weng et al.,, 2018)). Each ReLU layer can
be bounded as follows

$\displaystyle\underline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})\leq\,\sigma({\mathbf{x}}^{(i)})\,\leq\,\overline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})\,\
,$ (10)

$\forall\
{\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}$,
with
$\underline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})\coloneqq\underline{{\boldsymbol{\Sigma}}}^{(i)}{\mathbf{x}}^{(i)}+\underline{{\boldsymbol{\beta}}}^{(i)}$,
$\overline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})\coloneqq\overline{{\boldsymbol{\Sigma}}}^{(i)}{\mathbf{x}}^{(i)}+\overline{{\boldsymbol{\beta}}}^{(i)}$,
where the (diagonal) relaxation slope matrices
$\underline{{\boldsymbol{\Sigma}}}^{(i)},\overline{{\boldsymbol{\Sigma}}}^{(i)}$
and lower- and upper- offset vectors
$\underline{{\boldsymbol{\beta}}}^{(i)},\overline{{\boldsymbol{\beta}}}^{(i)}$
are given by

$\displaystyle\overline{{\boldsymbol{\Sigma}}}^{(i)}_{jj}\\!=\\!\begin{cases}\quad
0\\\ \quad 1\\\
\frac{{\boldsymbol{u}}^{(i)}_{j}}{{\boldsymbol{u}}^{(i)}_{j}-{\boldsymbol{\ell}}^{(i)}_{j}}\end{cases},\
\underline{{\boldsymbol{\Sigma}}}^{(i)}_{jj}\\!=\\!\begin{cases}\
0\\\\[2.84526pt] \ 1\\\\[2.84526pt] \,\alpha^{(i)}_{j}\end{cases},\
\overline{{\boldsymbol{\beta}}}^{(i)}_{j}=\begin{cases}\quad 0\\\ \quad 0\\\
\frac{-{\boldsymbol{u}}^{(i)}_{j}{\boldsymbol{\ell}}^{(i)}_{j}}{{\boldsymbol{u}}^{(i)}_{j}-{\boldsymbol{\ell}}^{(i)}_{j}}\end{cases},\
\underline{{\boldsymbol{\beta}}}^{(i)}_{j}=\begin{cases}\ 0\\\\[2.84526pt] \
0\\\\[2.84526pt] \ 0\\\\[2.84526pt] \end{cases},\
\begin{matrix}&\,\textnormal{if}\,\,j\in\mathcal{I}^{(i)}_{-}\\\\[2.84526pt]
&\,\textnormal{if}\,\,j\in\mathcal{I}^{(i)}_{+}\\\\[2.84526pt]
&\,\textnormal{if}\,\,j\in\mathcal{I}^{(i)}\end{matrix}$ (11)

with $0\leq\alpha^{(i)}_{j}\leq 1$, and where $\mathcal{I}^{(i)}_{-}$,
$\mathcal{I}^{(i)}_{+}$, and $\mathcal{I}^{(i)}$ denote the sets of
activations $j\in\llbracket n_{i}\rrbracket$ in layer $i$ where the lower and
upper pre-activation bounds
${\boldsymbol{\ell}}^{(i)}_{j},{\boldsymbol{u}}^{(i)}_{j}$ are both negative,
both positive, or span zero respectively. See Figure 2 (right) for an
illustration.

Proofs can be found in (Weng et al.,, 2018; Zhang et al.,, 2018; Salman et
al.,, 2019). See (Zhang et al.,, 2018) or (Liu et al.,, 2019) for how to
define the relaxation slope matrices, lower- and upper- offset vectors for
general activation functions.

Figure 2:  (Left) ReLU with pre-activation bounds
${\boldsymbol{\ell}}^{(i)}_{j}$ and ${\boldsymbol{u}}^{(i)}_{j}$. (Middle)
Optimal ReLU relaxation. (Right) Special relaxation where
$\underline{{\boldsymbol{\Sigma}}}^{(i)}=\overline{{\boldsymbol{\Sigma}}}^{(i)}$.
The slopes, upper- and lower- offsets are determined by
$\underline{{\boldsymbol{\Sigma}}}^{(i)}_{jj},\overline{{\boldsymbol{\Sigma}}}^{(i)}_{jj}$
and
$\underline{{\boldsymbol{\beta}}}^{(i)}_{j},\overline{{\boldsymbol{\beta}}}^{(i)}_{j}$
respectively.

Note that the linear upper bound is the optimal convex relaxation for the ReLU
non-linearity, cf. Equation 8. For the lower bound, the optimal convex
relaxation is not achievable as one linear function, however, and we can use
any “sub-gradient” relaxation
$\sigma_{j}({\mathbf{x}}^{(i)})=\alpha^{(i)}_{j}{\mathbf{x}}^{(i)}_{j}$ with
$0\leq\alpha^{(i)}_{j}\leq 1$, cf. Equation 8.

Choosing the same slope for the lower- and upper- relaxation functions, i.e.
$\underline{{\boldsymbol{\Sigma}}}^{(i)}=\overline{{\boldsymbol{\Sigma}}}^{(i)}$,
recovers Fast-Lin (Weng et al.,, 2018) and is equivalent to DeepZ (Singh et
al.,, 2018). The same procedure is also used to compute pre-activation bounds
in (Wong et al.,, 2018) (Algorithm 1). The case where the slopes
$\underline{{\boldsymbol{\Sigma}}}^{(i)}\neq\overline{{\boldsymbol{\Sigma}}}^{(i)}$
of the lower relaxation function
$\underline{{\boldsymbol{\Phi}}}^{(i)}(\cdot)$ are selected adaptively,
$\alpha^{(i)}_{j}\in\\{0,1\\}$, depending on which relaxation has the smaller
volume, recovers CROWN (Zhang et al.,, 2018) and is equivalent to DeepPoly
(Singh et al., 2019b, ).

#### Greedy Backsubstitution

The relaxation-based certification problem with linear bounding functions can
be solved greedily in a layer-by-layer fashion (Wong and Kolter,, 2018; Weng
et al.,, 2018; Zhang et al.,, 2018; Salman et al.,, 2019). For instance, to
obtain bounds on
${\mathbf{x}}^{(i+1)}_{j}={\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\sigma({\mathbf{x}}^{(i)})+{\mathbf{b}}^{(i+1)}_{j}$,
we greedily replace the non-linearities $\sigma({\mathbf{x}}^{(i)})$ with
their lower and upper relaxation functions,
$\underline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})$ resp.
$\overline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})$, in such a way that
we under-estimate the lower bounds ${\boldsymbol{\ell}}^{(i+1)}_{j}$ and over-
estimate the upper bounds ${\boldsymbol{u}}^{(i+1)}_{j}$. Specifically, the
bounding functions for $\sigma({\mathbf{x}}^{(i)})$ are chosen based on the
signs of the elements of the $j$-th row of the weight matrix,
${\mathbf{W}}^{(i+1)}_{j\cdot\cdot}$, i.e.

$\displaystyle{\boldsymbol{\ell}}^{(i+1)}_{j}$
$\displaystyle\geq\Big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\underline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\overline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})+{\mathbf{b}}^{(i)}_{j}$
(12)
$\displaystyle=\Big{(}\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\underline{{\boldsymbol{\Sigma}}}^{(i)}+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\overline{{\boldsymbol{\Sigma}}}^{(i)}\Big{)}\,{\mathbf{x}}^{(i)}$
$\displaystyle\quad+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\underline{{\boldsymbol{\beta}}}^{(i)}+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\overline{{\boldsymbol{\beta}}}^{(i)}+{\mathbf{b}}^{(i)}_{j}$
$\displaystyle{\boldsymbol{u}}^{(i+1)}_{j}$
$\displaystyle\leq\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\underline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\overline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})+{\mathbf{b}}^{(i)}_{j}$
$\displaystyle=\Big{(}\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\underline{{\boldsymbol{\Sigma}}}^{(i)}+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\overline{{\boldsymbol{\Sigma}}}^{(i)}\Big{)}\,{\mathbf{x}}^{(i)}$
$\displaystyle\quad+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\underline{{\boldsymbol{\beta}}}^{(i)}+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\overline{{\boldsymbol{\beta}}}^{(i)}+{\mathbf{b}}^{(i)}_{j}$

Similarly, in the expression for
${\mathbf{x}}^{(i)}={\mathbf{W}}^{(i)}\sigma({\mathbf{x}}^{(i-1)})+{\mathbf{b}}^{(i)}$
in the bounds above, we can greedily replace the non-linearities
$\sigma({\mathbf{x}}^{(i-1)})$ with their relaxation functions
$\underline{{\boldsymbol{\Phi}}}^{(i-1)}({\mathbf{x}}^{(i-1)})$ resp.
$\overline{{\boldsymbol{\Phi}}}^{(i-1)}({\mathbf{x}}^{(i-1)})$ depending on
the signs of
$\Big{(}\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\underline{{\boldsymbol{\Sigma}}}^{(i)}+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\overline{{\boldsymbol{\Sigma}}}^{(i)}\Big{)}{\mathbf{W}}^{(i)}$
and
$\Big{(}\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{-}\underline{{\boldsymbol{\Sigma}}}^{(i)}+\big{[}{\mathbf{W}}^{(i+1)}_{j\cdot\cdot}\big{]}_{+}\overline{{\boldsymbol{\Sigma}}}^{(i)}\Big{)}{\mathbf{W}}^{(i)}$,
thus obtaining linear bounds for ${\mathbf{x}}^{(i+1)}_{j}$ in terms of
${\mathbf{x}}^{(i-1)}$. By the same argument we can continue this
backsubstitution process until we reach the input ${\mathbf{x}}^{(0)}$, thus
getting bounds on ${\mathbf{x}}^{(i+1)}_{j}$ in terms of ${\mathbf{x}}^{(0)}$
of the following form, which holds for all
${\mathbf{x}}^{(0)}\in\mathcal{B}(\mathbf{{\mathbf{x}}})$,

$\displaystyle\underline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)}{\mathbf{x}}^{(0)}+\underline{{\boldsymbol{\gamma}}}_{j}^{(i+1)}\leq{\mathbf{x}}^{(i+1)}_{j}\leq\overline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)}{\mathbf{x}}^{(0)}+\overline{{\boldsymbol{\gamma}}}_{j}^{(i+1)}$
(13)

where
$\underline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)},\overline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)}$
capture the products of weight matrices and relaxation slope matrices, while
$\underline{{\boldsymbol{\gamma}}}_{j}^{(i+1)},\overline{{\boldsymbol{\gamma}}}_{j}^{(i+1)}$
collect products of weight matrices, relaxation matrices and bias terms.
Explicit expressions for
$\underline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)},\overline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)},\underline{{\boldsymbol{\gamma}}}_{j}^{(i+1)},\overline{{\boldsymbol{\gamma}}}_{j}^{(i+1)}$
can be found in (Weng et al.,, 2018; Zhang et al.,, 2018). The above recursion
is quite remarkable, as it allows to linearly bound the output of the non-
linear mapping
${\mathbf{x}}^{(i+1)}\equiv{\mathbf{x}}^{(i+1)}({\mathbf{x}}^{(0)})$ for all
${\mathbf{x}}^{(0)}\in\mathcal{B}({\mathbf{x}})$.

With the above expressions for the linear functions bounding
${\mathbf{x}}^{(i+1)}_{j}$ in terms of ${\mathbf{x}}^{(0)}$, we can compute
lower and upper bounds
${\boldsymbol{\ell}}^{(i+1)}_{j},{\boldsymbol{u}}^{(i+1)}_{j}$ by considering
the worst-case ${\mathbf{x}}^{(0)}\in\mathcal{B}({\mathbf{x}})$. For instance,
when the uncertainty set is an $\ell_{p}$-norm ball
$\mathcal{B}_{p}({\mathbf{x}};\epsilon)\coloneqq\\{{\mathbf{x}}^{*}:||{\mathbf{x}}^{*}-{\mathbf{x}}||_{p}\leq\epsilon\\}$,
the bounds are given as

$\displaystyle{\boldsymbol{\ell}}^{(i+1)}_{j}$
$\displaystyle\geq\underline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)}{\mathbf{x}}-\epsilon||\underline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)}||_{p^{*}}+\underline{{\boldsymbol{\gamma}}}_{j}^{(i+1)}$
(14) $\displaystyle{\boldsymbol{u}}^{(i+1)}_{j}$
$\displaystyle\leq\overline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)}{\mathbf{x}}+\epsilon||\overline{{\boldsymbol{\Lambda}}}_{j\cdot\cdot}^{(i+1)}||_{p^{*}}+\overline{{\boldsymbol{\gamma}}}_{j}^{(i+1)}$

where $p^{*}$ dentos the Hölder conjugate of p, given by $1/p+1/p^{*}=1$.

In practice, in order to be able to replace all the non-linearities with their
bounding functions from layer $i$ all the way down to the input layer, we need
pre-activation bounds for all layers $i^{\prime}<i$, since the expressions for
the relaxation slope matrices
$\underline{{\boldsymbol{\Sigma}}}^{(i^{\prime})},\overline{{\boldsymbol{\Sigma}}}^{(i^{\prime})}$
and offset vectors
$\underline{{\boldsymbol{\beta}}}^{(i^{\prime})},\overline{{\boldsymbol{\beta}}}^{(i^{\prime})}$
depend on those pre-activation bounds. The full greedy solution thus proceeds
in a layer-by-layer fashion, starting from the first layer up to the last
layer, where for each layer the backsubstitution to the input is computed
based on the pre-activation bounds of previous layers (computed with the same
greedy approach). Hence, the computational complexity of the full greedy
solution of the relaxation-based certification problem is quadratic in the
number of layers $L$ of the network. The complete algorithm is shown in
Section 5.3 in the Appendix.

Finally, note that while the activation over-approximation introduces
looseness, relaxation-based bounds admit cancellations between positive and
negative entries in weight matrices that are otherwise missed when considering
the positive and negative parts of the weight matrices separately as in the
Interval Bound Propagation algorithm.

#### Dual formulations

Instead of solving the certification problem
$\mathcal{O}(\mathbf{c},c_{0},L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$ in
Equation 5 or the corresponding relaxation
$\mathcal{C}(\mathbf{c},c_{0},L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$ in
Equation 6 in their primal forms, we can also solve the corresponding dual
formulations. The Lagrangian dual of the certification problem in Equation 5
is given by (Dvijotham et al.,, 2018; Salman et al.,, 2019)

$\displaystyle g_{\mathcal{O}}({\boldsymbol{\mu}}^{\llbracket
L\rrbracket})=\min_{{\mathbf{x}}^{\llbracket 0,L\rrbracket}}$ $\displaystyle\
\mathbf{c}^{\top}{\mathbf{x}}^{(L)}+c_{0}+\sum_{i=1}^{L}{\boldsymbol{\mu}}^{(i)\top}\left({\mathbf{x}}^{(i)}-{\mathbf{W}}^{(i)}\sigma({\mathbf{x}}^{(i-1)})-{\mathbf{b}}^{(i)}\right)$
(15) s.t. $\displaystyle\ {\mathbf{x}}^{(0)}\in\mathcal{B}({\mathbf{x}})\,,\
{\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket L\\!-\\!1\rrbracket$

The Lagrangian dual of the relaxation-based certification problem in Equation
6 is given by (Wong and Kolter,, 2018; Salman et al.,, 2019)

$\displaystyle g_{\mathcal{C}}({\boldsymbol{\mu}}^{\llbracket
L\rrbracket},\underline{\boldsymbol{\lambda}}^{\llbracket
0,L\\!-\\!1\rrbracket},\overline{\boldsymbol{\lambda}}^{\llbracket
0,L\\!-\\!1\rrbracket})$ (16) $\displaystyle=\min_{{\mathbf{x}}^{\llbracket
0,L\rrbracket},{\mathbf{z}}^{\llbracket 0,L\\!-\\!1\rrbracket}}$
$\displaystyle\mathbf{c}^{\top}{\mathbf{x}}^{(L)}+c_{0}+\sum_{i=1}^{L}{\boldsymbol{\mu}}^{(i)\top}\\!\left({\mathbf{x}}^{(i)}-{\mathbf{W}}^{(i)}{\mathbf{z}}^{(i-1)}-{\mathbf{b}}^{(i)}\right)$
$\displaystyle-\sum_{i=0}^{L-1}\underline{\boldsymbol{\lambda}}^{(i)\top}\left({\mathbf{z}}^{(i)}-\underline{\sigma}({\mathbf{x}}^{(i)})\right)+\sum_{i=0}^{L-1}\overline{\boldsymbol{\lambda}}^{(i)\top}\left({\mathbf{z}}^{(i)}-\overline{\sigma}({\mathbf{x}}^{(i)})\right)$
s.t. $\displaystyle\ {\mathbf{x}}^{(0)}\in\mathcal{B}({\mathbf{x}})\,,\
{\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket L\\!-\\!1\rrbracket$

where ${\mathbf{z}}^{(i)}$ represent the post-activation variables.

By weak duality (Boyd et al.,, 2004), we have for the original dual

$\displaystyle
d^{*}_{\mathcal{O}}\coloneqq\max_{{\boldsymbol{\mu}}^{\llbracket
L\rrbracket}}g_{\mathcal{O}}({\boldsymbol{\mu}}^{\llbracket L\rrbracket})\leq
p^{*}_{\mathcal{O}}$ (17)

respectively for the convex-relaxation based dual

$\displaystyle
d^{*}_{\mathcal{C}}\coloneqq\max_{{\boldsymbol{\mu}}^{\llbracket
L\rrbracket},\underline{\boldsymbol{\lambda}}^{\llbracket
0,L\\!-\\!1\rrbracket}\geq 0,\overline{\boldsymbol{\lambda}}^{\llbracket
0,L\\!-\\!1\rrbracket}\geq 0}g_{\mathcal{C}}({\boldsymbol{\mu}}^{\llbracket
L\rrbracket},\underline{\boldsymbol{\lambda}}^{\llbracket
0,L\\!-\\!1\rrbracket},\overline{\boldsymbol{\lambda}}^{\llbracket
0,L\\!-\\!1\rrbracket})\leq p^{*}_{\mathcal{C}}$ (18)

Hence, if $d^{*}_{\mathcal{O}}>0$ resp. $d^{*}_{\mathcal{C}}>0$ for all
$\mathbf{c}\in\\{\mathbf{e}_{\hat{k}({\mathbf{x}})}-\mathbf{e}_{k}\\}_{k\neq\hat{k}({\mathbf{x}})}$,
with $c_{0}=0$ and valid pre-activation bounds
${\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket}$, the
network is certifiably robust with respect to ${\mathbf{x}}$ and
$\mathcal{B}({\mathbf{x}})$.

In fact, for the convex-relaxation based certification problem, one can show
that strong duality ($d^{*}_{\mathcal{C}}=p^{*}_{\mathcal{C}}$) holds under
relatively mild conditions (finite Lipschitz constant for the bounding
functions $\underline{\sigma}(\cdot)$, $\overline{\sigma}(\cdot)$), see
Theorem 4.1 in (Salman et al.,, 2019). Moreover, one can even show that the
dual of the optimal single-neuron convex relaxation based certification
problem is equivalent to the dual of the original certification problem, i.e.
$d^{*}_{\mathcal{C}_{\textnormal{opt}}}=d^{*}_{\mathcal{O}}$, see Theorem 4.2
in (Salman et al.,, 2019).

Despite these equivalences, there are still good reasons to solve the dual
instead of the primal problem. Salman et al., (2019) recommend solving the
dual problem because (i) the dual problem can be formulated as an
unconstrained optimization problem, whereas the primal is a constrained
optimization problem and (ii) the dual optimization process can be stopped
anytime to give a valid lower bound on $p^{*}_{\mathcal{O}}$ (thanks to weak
duality).

#### Single-Neuron Relaxation Barrier

The effectiveness of existing single-neuron relaxation-based verifiers is
inherently limited by the tightness of the optimal univariate relaxation
(Salman et al.,, 2019). In a series of highly compute-intense experiments,
Salman et al., (2019) found that optimal single-neuron relaxation based
verification, i.e. solving the certification problems
$\mathcal{C}_{\textrm{opt}}(\pm\mathbf{e}_{j},0,i,{\boldsymbol{\ell}}^{\llbracket
i\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket i\\!-\\!1\rrbracket})$ for
all neurons $j\in\llbracket n_{i}\rrbracket$ and layers $i$, does not
significantly improve upon the gap between verifiers that greedily compute
approximate pre-activation bounds and only solve the relaxed (primal or dual)
certification problem
$\mathcal{C}_{\textrm{opt}}(\mathbf{c},0,L,{\boldsymbol{\ell}}^{\llbracket
L\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket L\\!-\\!1\rrbracket})$ with
$\mathbf{c}\in\\{\mathbf{e}_{\hat{k}({\mathbf{x}})}-\mathbf{e}_{k}\\}_{k\neq\hat{k}({\mathbf{x}})}$
for the logit layer, and the exact mixed integer linear programming (MILP)
verifier from (Tjeng et al.,, 2018), suggesting that there is an inherent
barrier to tight verification for single neuron relaxations. To improve the
tightness of relaxation-based certification methods we therefore have to go
beyond single-neuron relaxations.

##  3 Multi-Neuron Relaxation

The single-neuron relaxation barrier can be bypassed by considering multi-
neuron relaxations, i.e. by approximating multiple non-linearities jointly
instead of separately, as suggested in the recent LP-solver based “k-ReLU”
relaxation framework of Singh et al., 2019a . Singh et al.’s multi-neuron
relaxation framework is specific to networks with ReLU non-linearities but can
otherwise be incorporated into any certification method (primal or dual) that
operates with pre-activation bounds, including (Weng et al.,, 2018; Zhang et
al.,, 2018; Wong and Kolter,, 2018).

In this Section, we show how the “k-ReLU” relaxation framework obtains tighter
activation bounds by leveraging additional relational constraints among groups
of neurons. In particular, we show how additional pre-activation bounds can be
mapped to corresponding post-activation bounds and how they can in turn be
used to obtain tighter bounds in higher layers, as illustrated in Figure 3. By
capturing interactions between neurons, the k-ReLU method is able to overcome
the convex barrier imposed by the single neuron relaxation (Salman et al.,,
2019).

Figure 3:  Multi-neuron relaxations obtain tighter correlation-aware
activation bounds by leveraging additional relational constraints among groups
of neurons, thereby overcoming the convex barrier imposed by the single-neuron
relaxations.

Singh et al.’s experimental results indicate that k-ReLU enables significantly
more precise certification than existing state-of-the-art verifiers while
maintaining scalability. To illustrate the precision gain, Singh et al.
measure the volume of the output bounding box computed after propagating an
$\ell_{\infty}$-ball of radius $\epsilon=0.015$ through a fully connected
network with $9$ layers containing $200$ neurons each. They find that the
volume of the output from 2-ReLU resp. 3-ReLU relaxation is 7 resp. 9 orders
of magnitude smaller than from single-neuron relaxation-based DeepPoly
verifier (Singh et al., 2019a, ). Similarly, on the $9\times 200$ fully-
connected network resp. a convolutional network, k-ReLU certifies $506$ resp.
$347$ adversarial regions whereas the single-neuron relaxation based
RefineZono verifier certifies $316$ resp. $179$ adversarial regions (Singh et
al., 2019a, ).

|  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 4: Visualization of different relaxation methods. (a) Uncertainty set
after the first affine layer (i.e. $\mathcal{S}^{(1)}$ projected onto
${\mathbf{x}}^{(1)}$) and its ReLU projection, (b) uncertainty set, its ReLU
projection, single-neuron relaxation (SNR) with same relaxation slope matrices
$\underline{{\boldsymbol{\Sigma}}}^{(1)}=\overline{{\boldsymbol{\Sigma}}}^{(1)}$
and SNR with optimal single-neuron relaxation (note that the same slope SNR
set contains the optimal SNR set which in turn contains the ReLU projection
set), (c) uncertainty set, its ReLU projection and the k-ReLU relaxation
(corresponding to the convex hull of the ReLU projection), (d) everything
combined for comparison. Note that corresponding pre- and post-activation
pairs are superimposed on the same axis, i.e. both ${\mathbf{x}}^{(1)}_{1}$
and ${\mathbf{z}}^{(1)}_{1}$ are plotted (superimposed) on the horizontal axis
while both ${\mathbf{x}}^{(1)}_{2}$ and ${\mathbf{z}}^{(1)}_{2}$ are plotted
(superimposed) on the vertical axis. This style of visialization, where pre-
and post-activation pairs are superimposed, allows to easily convey the effect
of the non-linearities.

#### ReLU branch polytopes

Following Singh et al., 2019a , we consider the pre-activations
${\mathbf{x}}^{(i)}$ and the post-activations ${\mathbf{z}}^{(i)}$ as separate
neurons. Let $\mathcal{S}^{(\ell)}$ be a convex set computed via some
relaxation based certification method approximating the set of values that
neurons $({\mathbf{x}}^{\llbracket 0,\ell\rrbracket},{\mathbf{z}}^{\llbracket
0,\ell-1\rrbracket})$, including the pre-activations ${\mathbf{x}}^{(\ell)}$
but excluding the post-activations ${\mathbf{z}}^{(\ell)}$, can take with
respect to $\mathcal{B}({\mathbf{x}})$,

$\displaystyle\mathcal{S}^{(\ell)}\coloneqq\Big{\\{}$
$\displaystyle({\mathbf{x}}^{\llbracket
0,L\rrbracket},{\mathbf{z}}^{\llbracket 0,L\\!-1\\!\rrbracket})\ :\ $ (19)
$\displaystyle{\mathbf{x}}^{(i)}={\mathbf{W}}^{(i)}{\mathbf{z}}^{(i-1)}+{\mathbf{b}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket\ell\rrbracket$
$\displaystyle\underline{\sigma}({\mathbf{x}}^{(i)})\leq{\mathbf{z}}^{(i)}\leq\overline{\sigma}({\mathbf{x}}^{(i)})\,\
\textnormal{for}\ i\in\llbracket\ell\\!-\\!1\rrbracket$
$\displaystyle{\mathbf{x}}^{(0)}\in\mathcal{B}({\mathbf{x}})\,,\
{\boldsymbol{\ell}}^{(i)}\leq{\mathbf{x}}^{(i)}\leq{\boldsymbol{u}}^{(i)}\,\
\textnormal{for}\ i\in\llbracket\ell\\!-\\!1\rrbracket\Big{\\}}$

Note that, variables
$({\mathbf{x}}^{\llbracket\ell+1,L\rrbracket},{\mathbf{z}}^{\llbracket\ell,L\\!-\\!1\rrbracket})$
that don’t appear in the constraints, are considered to be unconstrained, i.e.
they take values in the entire real number line.

In general, pre-activations ${\mathbf{x}}^{(\ell)}_{j}$ can take both positive
and negative values in $\mathcal{S}^{(\ell)}$. For each ReLU activation where
the input, i.e. the corresponding pre-activation, can take both positive and
negative values, two branches have to be considered. Define the convex
polytopes induced by the two branches of the $j$-th ReLU unit at layer $\ell$
as

$\displaystyle C^{(\ell)}_{j\,+}$
$\displaystyle\coloneqq\\{({\mathbf{x}}^{\llbracket
0,L\rrbracket},{\mathbf{z}}^{\llbracket
0,L\\!-1\\!\rrbracket}):{\mathbf{x}}^{(\ell)}_{j}\geq 0,\
{\mathbf{z}}^{(\ell)}_{j}={\mathbf{x}}^{(\ell)}_{j}\\}$ (20) $\displaystyle
C^{(\ell)}_{j\,-}$ $\displaystyle\coloneqq\\{({\mathbf{x}}^{\llbracket
0,L\rrbracket},{\mathbf{z}}^{\llbracket
0,L\\!-1\\!\rrbracket}):{\mathbf{x}}^{(\ell)}_{j}\leq 0,\
{\mathbf{z}}^{(\ell)}_{j}=0\\}\,,$

which can be written more concisely as

$\displaystyle C^{(\ell)}_{j\,s_{j}}$
$\displaystyle\coloneqq\Big{\\{}({\mathbf{x}}^{\llbracket
0,L\rrbracket},{\mathbf{z}}^{\llbracket
0,L\\!-1\\!\rrbracket}):s_{j}{\mathbf{x}}^{(\ell)}_{j}\geq 0,\
{\mathbf{z}}^{(\ell)}_{j}=\frac{1+s_{j}}{2}{\mathbf{x}}^{(\ell)}_{j}\Big{\\}}\,.$
(21)

Next, we introduce some notation to describe the polytopes representing the
possible ways of selecting one ReLU branch per neuron. Let
$J\subseteq\llbracket n_{\ell}\rrbracket$ be some index set over neurons, with
cardinality $|J|$. For a specific configuration
$(s_{j_{1}},\dots,s_{j_{|J|}})\in\\{+,-\\}^{|J|}$ of individual ReLU branches,
out of all $2^{|J|}$ possible configurations, let

$\displaystyle Q^{(\ell)}_{J,(s_{j_{1}},\dots,s_{j_{|J|}})}$
$\displaystyle\coloneqq\bigcap_{j\in J}C^{(\ell)}_{j\,s_{j}}$ (22)

be the convex polytope defined as the intersection of the $|J|$ individual
ReLU branch polytopes $\\{C^{(\ell)}_{j\,s_{j}}\\}_{j\in J}$. Based on the
linear inequality description of the individual ReLU branch polytopes, this
can also be written explicitly as

$\displaystyle
Q^{(\ell)}_{J,(s_{j_{1}},\dots,s_{j_{|J|}})}=\Big{\\{}({\mathbf{x}}^{\llbracket
0,L\rrbracket},{\mathbf{z}}^{\llbracket
0,L\\!-1\\!\rrbracket}):s_{j}{\mathbf{x}}^{(\ell)}_{j}\geq 0,\
{\mathbf{z}}^{(\ell)}_{j}=\frac{1\\!+\\!s_{j}}{2}{\mathbf{x}}^{(\ell)}_{j}\,\
\textnormal{for}\ j\in J\Big{\\}}\,.$ (23)

Next, define $Q^{(\ell)}_{J}$ as the collection of all $2^{|J|}$ convex
polytopes $Q^{(\ell)}_{J,(s_{j_{1}},\dots,s_{j_{|J|}})}$, indexed by all
possible ReLU branch configurations
$(s_{j_{1}},\dots,s_{j_{|J|}})\in\\{+,-\\}^{|J|}$,

$\displaystyle Q^{(\ell)}_{J}\coloneqq\bigg{\\{}\bigcap_{j\in
J}C^{(\ell)}_{j\,s_{j}}\ \Big{|}\
(s_{j_{1}},\dots,s_{j_{|J|}})\in\\{+,-\\}^{|J|}\bigg{\\}}\,.$ (24)

#### Additional Relational Constraints

The k-ReLU relaxation framework bypasses the single-neuron convex barrier by
leveraging additional relational constraints among groups of neurons.
Specifically, the k-ReLU framework computes bounds on additional relational
constraints of the form $\sum_{j\in J}a_{j}{\mathbf{x}}^{(\ell)}_{j}$. These
additional relational constraints, together with the usual interval bounds,
are captured by the convex polytope $P^{(\ell)}_{J,\mathcal{A}}$. Formally,
let $P^{(\ell)}_{J,\mathcal{A}}\supseteq\mathcal{S}^{(\ell)}$ be a convex
polytope containing interval constraints (for $(a_{j_{1}},\dots,a_{j_{|J|}})$
with $a_{j}\in\\{-1,0,1\\}$, $\sum_{j\in J}|a_{j}|=1$) and relational
constraints (for general $(a_{j_{1}},\dots,a_{j_{|J|}})$, with
$a_{j}\in\mathbb{R}$) over neurons ${\mathbf{x}}^{(\ell)}_{j},j\in J$, defined
as

$\displaystyle
P^{(\ell)}_{J,\mathcal{A}}\coloneqq\Big{\\{}({\mathbf{x}}^{\llbracket
0,L\rrbracket},{\mathbf{z}}^{\llbracket 0,L\\!-1\\!\rrbracket}):$
$\displaystyle\sum_{j\in J}a_{j}{\mathbf{x}}^{(\ell)}_{j}\leq
c^{(\ell)}(a_{j_{1}},\dots,a_{j_{|J|}})\Big{|}(a_{j_{1}},\dots,a_{j_{|J|}})\in\mathcal{A}\Big{\\}}\,,$
(25)

where the set $\mathcal{A}$ contains the coefficient-tuples of all the
constraints defining $P^{(\ell)}_{J,\mathcal{A}}$.

In practice, the corresponding bounds
$c^{(\ell)}(a_{j_{1}},\dots,a_{j_{|J|}})$ on interval and relational
constraints over neurons ${\mathbf{x}}^{(\ell)}_{j},j\in J$, can be computed
using any of the existing bound computation algorithms, e.g. (Zhang et al.,,
2018), Algorithm 1 in (Wong and Kolter,, 2018) or DeepPoly (Singh et al.,
2019a, ). Using the notation for the relaxation-based certification problem in
Equation 6, the bounds are given as

$\displaystyle c^{(\ell)}(a_{j_{1}},\dots,a_{j_{|J|}})$
$\displaystyle=\mathcal{C}(-\mathbf{a},0,\ell,{\boldsymbol{\ell}}^{\llbracket\ell\\!-\\!1\rrbracket},{\boldsymbol{u}}^{\llbracket\ell\\!-\\!1\rrbracket})$
(26)

with $\mathbf{a}=\sum_{j\in J}a_{j}\mathbf{e}_{j}$, where $\mathbf{e}_{j}$
denotes the $j$-th canonical basis vector.

Ideally, one would like $P^{(\ell)}_{J,\mathcal{A}}$ to be the projection of
$\mathcal{S}^{(\ell)}$ onto the variables ${\mathbf{x}}^{(\ell)}_{j}$ indexed
by $J$. However, computing this projection is prohibitively expensive. Singh
et al., 2019a  heuristically found
$\mathcal{A}=\\{(a_{j_{1}},\dots,a_{j_{|J|}})\in\\{-1,0,1\\}^{|J|}\setminus(0,\dots,0)\\}$,
containing $3^{|J|}-1$ constraints ($2|J|$ interval and $3^{|J|}-2|J|-1$
relational), to work well in practice. It remains an open problem of whether
there exists a theoretically optimal arrangement of a given number of
additional relational constraints. Geometrically, $P^{(\ell)}_{J,\mathcal{A}}$
is an over-approximation of the projection of $\mathcal{S}^{(\ell)}$ onto the
variables ${\mathbf{x}}^{(\ell)}_{j}$ indexed by $J$.

#### Optimal Convex Multi-Neuron Relaxation

The optimal convex relaxation of the $n_{\ell}$ ReLU assignments considers all
$n_{\ell}$ neurons jointly

$\displaystyle\mathcal{S}^{(\ell)}_{\textnormal{opt}}$
$\displaystyle=\mathcal{CH}\Big{(}\bigcup_{Q\in Q^{(\ell)}_{\llbracket
n_{\ell}\rrbracket}}\\{\mathcal{S}^{(\ell)}\cap Q\\}\,\Big{)}\,,$ (27)

where $Q^{(\ell)}_{\llbracket n_{\ell}\rrbracket}$ is the collection of the
$2^{n_{\ell}}$ convex polytopes
$Q^{(\ell)}_{(s_{1},\dots,s_{n_{\ell}})}\coloneqq\bigcap_{j\in\llbracket
n_{\ell}\rrbracket}C^{(\ell)}_{j\,s_{j}}$ indexed by all possible ReLU branch
configurations across neurons at layer $\ell$, and where $\mathcal{CH}$
denotes the Convex-Hull. Note that $\mathcal{S}^{(\ell)}_{\textnormal{opt}}$
can equivalently be written as

$\displaystyle\mathcal{S}^{(\ell)}_{\textnormal{opt}}$
$\displaystyle=\mathcal{CH}\Big{(}\bigcup_{(s_{1},\cdots,s_{n_{\ell}})\in\\{+,-\\}^{n_{\ell}}}\\{\mathcal{S}^{(\ell)}\cap
Q^{(\ell)}_{(s_{1},\dots,s_{n_{\ell}})}\\}\,\Big{)}\,.$ (28)

Computing $\mathcal{S}^{(\ell)}_{\textnormal{opt}}$ is practically infeasible,
however, as the convex hull computation has exponential cost in the number of
neurons (Singh et al., 2019a, ).

#### k-ReLU Relaxation

The k-ReLU framework partitions the $n_{\ell}$ neurons into $n_{\ell}/k$
disjoint sub-groups of size $k$, then jointly relaxes the neurons within each
sub-group, as a compromise between the practically infeasible joint relaxation
of all $n_{\ell}$ neurons and the single-neuron relaxation barrier arising
when relaxing each neuron individually.

Suppose that $n_{\ell}$ is divisible by $k$. Let
$\mathcal{J}=\\{J_{i}\\}_{i=1}^{n_{\ell}/k}$ be a partition of the set of
neurons $\llbracket n_{\ell}\rrbracket$ such that each group
$J_{i}\in\mathcal{J}$ contains exactly $|J_{i}|=k$ indices. Singh et al.,
2019a  k-ReLU framework computes the following convex relaxation

$\displaystyle\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}$
$\displaystyle=\mathcal{S}^{(\ell)}\cap\bigcap_{i=1}^{n_{\ell}/k}\mathcal{CH}\Big{(}\bigcup_{Q\in
Q^{(\ell)}_{J_{i}}}\\{P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}\cap
Q\\}\,\Big{)}\,.$ (29)

Conceptually, each convex polytope $P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}$
(containing interval and relational constraints on neurons
${\mathbf{x}}^{(\ell)}_{j},j\in J_{i}$) is intersected with convex ReLU branch
polytopes $Q$ from the set $Q^{(\ell)}_{J_{i}}$, producing $2^{k}$ convex
polytopes $\\{P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}\cap Q\\}_{{Q\in
Q^{(\ell)}_{J_{i}}}}$, one for each possible branch ${Q\in
Q^{(\ell)}_{J_{i}}}$.

The union $\bigcup_{Q\in
Q^{(\ell)}_{J_{i}}}\\{P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}\cap Q\\}$ of all
$2^{k}$ convex polytopes captures222“over-approximates” (to be more precise),
since $P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}$ is an over-approximation of the
projection of $\mathcal{S}^{(\ell)}$ onto the variables indexed by $J_{i}$.
the uncertainty set over post-activations ${\mathbf{z}}^{(\ell)}_{j}$ indexed
by $j\in J_{i}$. As the different groups of neurons are disjoint $J_{i}\cap
J_{j}=\emptyset$, we can combine those group-specific post-activation
uncertainty sets by intersecting their convex hulls
$\mathcal{K}_{i}=\mathcal{CH}\big{(}\bigcup_{Q\in
Q^{(\ell)}_{J_{i}}}\\{P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}\cap Q\\}\,\big{)}$.
From this, $\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}$ is obtained by
intersection with the convex set $\mathcal{S}^{(\ell)}$. An illustration of
the k-ReLU relaxation can be found in Figure 4.

Singh et al., 2019a  heuristically chose the partition
$\mathcal{J}=\\{J_{i}\\}_{i=1}^{n_{\ell}/k}$ such that neurons $j\in\llbracket
n_{\ell}\rrbracket$ are grouped according to the area of their triangle
relaxation (i.e. neurons with similar triangle relaxation areas are grouped
together). It remains an open problem whether there exists a theoretically
optimal partitioning of the groups of neurons.

Finally, we note that $\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}$ can
equivalently be written as

$\displaystyle\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}$
$\displaystyle=\mathcal{S}^{(\ell)}\cap\bigcap_{i=1}^{n_{\ell}/k}\mathcal{CH}\Big{(}\bigcup_{\begin{subarray}{c}(s_{j_{1}},\dots,s_{j_{k}})\in\\{+,-\\}^{|J_{i}|}\\\
j_{1},\dots,j_{k}\in
J_{i}\end{subarray}}\\{P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}\cap
Q^{(\ell)}_{\begin{subarray}{c}J_{i},\,(s_{j_{1}},\dots,s_{j_{k}})\end{subarray}}\\}\,\Big{)}\,.$
(30)

#### Convex Hull Computation

Singh et al. use the cdd library to compute convex hulls (cdd,, 2021). cdd is
an implementation of the double description method by Motzkin, that allows to
compute all vertices (i.e. extreme points) and extreme rays of a general
convex polytope given as a system of linear inequalities. cdd also implements
the reverse operation, allowing to compute the convex hull from a set of
vertices. In practice, cdd is used first to compute the vertices of the convex
polytopes $\\{P^{(\ell)}_{J_{i},\,\mathcal{A}_{i}}\cap Q\\}_{{Q\in
Q^{(\ell)}_{J_{i}}}}$ and second to compute the convex hull of the union of
the vertices of all these polytopes.

#### Tighter bounds in higher layers

Finally, we show how the k-ReLU relaxation
$\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}$ can be used to obtain tighter
bounds in higher layers. The k-ReLU method computes _refined_ pre-activation
bounds ${\boldsymbol{\ell}}^{(i+1)}_{j},{\boldsymbol{u}}^{(i+1)}_{j}$ for
neurons at layer $\ell+1$, by maximizing and minimizing
${\mathbf{x}}^{(i+1)}_{j}={\mathbf{W}}^{(\ell+1)}_{j\cdot\cdot}{\mathbf{z}}^{(\ell)}+{\mathbf{b}}^{(\ell+1)}_{j}$
with respect to ${\mathbf{z}}^{(\ell)}$ subject to
${\mathbf{z}}^{(\ell)}\in\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}$, i.e.

$\displaystyle{\boldsymbol{\ell}}^{(\ell+1)}_{j}=\min_{{\mathbf{z}}^{(\ell)}\in\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}}{\mathbf{W}}^{(\ell+1)}_{j\cdot\cdot}{\mathbf{z}}^{(\ell)}+{\mathbf{b}}^{(\ell+1)}_{j}$
(31)
$\displaystyle{\boldsymbol{u}}^{(\ell+1)}_{j}=\max_{{\mathbf{z}}^{(\ell)}\in\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}}{\mathbf{W}}^{(\ell+1)}_{j\cdot\cdot}{\mathbf{z}}^{(\ell)}+{\mathbf{b}}^{(\ell+1)}_{j}$

Since all the constraints in $\mathcal{S}^{(\ell)}_{\textnormal{k-ReLU}}$ are
linear, we can use an LP-solver for the maximization and minimization. Singh
et al. use the gurobi solver (Gurobi,, 2021).

#### Interpretation

Computing bounds on the additional relational constraints corresponds to
solving a bounding problem on a widened network, that is equal to the original
network up to layer $\ell-1$ but with a “modified” $\ell$-th layer which
includes additional neurons representing linear combinations of rows of the
$\ell$-th layer weight matrix with coefficients determined by
$(a_{j_{1}},\dots,a_{j_{|J|}})$. See Figure 3 for an illustration.

##  4 Discussion

We have shown how the recent “k-ReLU” framework of Singh et al., 2019a
obtains tighter correlation-aware activation bounds by leveraging additional
relational constraints among groups of neurons. In particular, we have shown
how additional pre-activation constraints can be mapped to corresponding post-
activation constraints and how they can in turn be used to obtain tighter pre-
activation bounds in higher layers. The k-ReLU framework is specific to ReLU
networks but can otherwise be incorporated into any certification method that
operates with pre-activation bounds.

The main degrees of freedom in k-ReLU are the partitioning of the neurons into
sub-groups and the choice of coefficients $a_{j}$ in the additional relational
constraints. We consider it to be an interesting future avenue of research to
investigate whether there is a theoretically optimal partitioning as well as
choice for the coefficients of the additiona relational constraints.

## Acknowledgments

I would like to thank Gagandeep Singh for the insightful discussions we had
and the helpful comments on the k-ReLU method he provided.

## References

  * Athalye et al., (2018) Athalye, A., Carlini, N., and Wagner, D. (2018).  Obfuscated gradients give a false sense of security: Circumventing defenses to adversarial examples.  In International Conference on Machine Learning, pages 274–283.
  * Biggio et al., (2013) Biggio, B., Corona, I., Maiorca, D., Nelson, B., Srndic, N., Laskov, P., Giacinto, G., and Roli, F. (2013).  Evasion attacks against machine learning at test time.  In Joint European conference on machine learning and knowledge discovery in databases, pages 387–402. Springer.
  * Boyd et al., (2004) Boyd, S., Boyd, S. P., and Vandenberghe, L. (2004).  Convex optimization.  Cambridge university press.
  * Carlini et al., (2017) Carlini, N., Katz, G., Barrett, C., and Dill, D. L. (2017).  Provably minimally-distorted adversarial examples.  arXiv preprint arXiv:1709.10207.
  * (5) Carlini, N. and Wagner, D. (2017a).  Adversarial examples are not easily detected: Bypassing ten detection methods.  In Proceedings of the 10th ACM Workshop on Artificial Intelligence and Security, pages 3–14. ACM.
  * (6) Carlini, N. and Wagner, D. (2017b).  Towards evaluating the robustness of neural networks.  In 2017 IEEE Symposium on Security and Privacy (SP), pages 39–57. IEEE.
  * cdd, (2021) cdd (2021).  pycddlib, https://pypi.org/project/pycddlib/.
  * Cheng et al., (2017) Cheng, C.-H., Nührenberg, G., and Ruess, H. (2017).  Maximum resilience of artificial neural networks.  In International Symposium on Automated Technology for Verification and Analysis, pages 251–268. Springer.
  * Dvijotham et al., (2018) Dvijotham, K., Stanforth, R., Gowal, S., Mann, T. A., and Kohli, P. (2018).  A dual approach to scalable verification of deep networks.  In UAI, volume 1, page 2.
  * Ehlers, (2017) Ehlers, R. (2017).  Formal verification of piece-wise linear feed-forward neural networks.  In International Symposium on Automated Technology for Verification and Analysis, pages 269–286. Springer.
  * Feinman et al., (2017) Feinman, R., Curtin, R. R., Shintre, S., and Gardner, A. B. (2017).  Detecting adversarial samples from artifacts.  arXiv preprint arXiv:1703.00410.
  * Goodfellow et al., (2014) Goodfellow, I. J., Shlens, J., and Szegedy, C. (2014).  Explaining and harnessing adversarial examples.  arXiv preprint arXiv:1412.6572.
  * Grosse et al., (2017) Grosse, K., Manoharan, P., Papernot, N., Backes, M., and McDaniel, P. (2017).  On the (statistical) detection of adversarial examples.  arXiv preprint arXiv:1702.06280.
  * Gurobi, (2021) Gurobi (2021).  Gurobi optimizer, http://www.gurobi.com.
  * Hein and Andriushchenko, (2017) Hein, M. and Andriushchenko, M. (2017).  Formal guarantees on the robustness of a classifier against adversarial manipulation.  In Advances in Neural Information Processing Systems, pages 2266–2276.
  * Katz et al., (2017) Katz, G., Barrett, C., Dill, D. L., Julian, K., and Kochenderfer, M. J. (2017).  Reluplex: An efficient smt solver for verifying deep neural networks.  In International Conference on Computer Aided Verification, pages 97–117. Springer.
  * Kurakin et al., (2016) Kurakin, A., Goodfellow, I., and Bengio, S. (2016).  Adversarial examples in the physical world.  arXiv preprint arXiv:1607.02533.
  * Liu et al., (2019) Liu, C., Tomioka, R., and Cevher, V. (2019).  On certifying non-uniform bound against adversarial attacks.  arXiv preprint arXiv:1903.06603.
  * Lomuscio and Maganti, (2017) Lomuscio, A. and Maganti, L. (2017).  An approach to reachability analysis for feed-forward relu neural networks.  arXiv preprint arXiv:1706.07351.
  * Madry et al., (2018) Madry, A., Makelov, A., Schmidt, L., Tsipras, D., and Vladu, A. (2018).  Towards deep learning models resistant to adversarial attacks.  In International Conference on Learning Representations.
  * Metzen et al., (2017) Metzen, J. H., Genewein, T., Fischer, V., and Bischoff, B. (2017).  On detecting adversarial perturbations.  arXiv preprint arXiv:1702.04267.
  * Mirman et al., (2018) Mirman, M., Gehr, T., and Vechev, M. (2018).  Differentiable abstract interpretation for provably robust neural networks.  In International Conference on Machine Learning, pages 3578–3586.
  * Müller et al., (2021) Müller, M. N., Makarchuk, G., Singh, G., Püschel, M., and Vechev, M. (2021).  Prima: Precise and general neural network certification via multi-neuron convex relaxations.  arXiv preprint arXiv:2103.03638.
  * Qin et al., (2019) Qin, C., O’Donoghue, B., Bunel, R., Stanforth, R., Gowal, S., Uesato, J., Swirszcz, G., Kohli, P., et al. (2019).  Verification of non-linear specifications for neural networks.  arXiv preprint arXiv:1902.09592.
  * Raghunathan et al., (2018) Raghunathan, A., Steinhardt, J., and Liang, P. (2018).  Certified defenses against adversarial examples.  arXiv preprint arXiv:1801.09344.
  * Roth et al., (2019) Roth, K., Kilcher, Y., and Hofmann, T. (2019).  The odds are odd: A statistical test for detecting adversarial examples.  In International Conference on Machine Learning, pages 5498–5507.
  * Roth et al., (2020) Roth, K., Kilcher, Y., and Hofmann, T. (2020).  Adversarial training is a form of data-dependent operator norm regularization.  In Advances in Neural Information Processing Systems 34, pages 2015–2025.
  * Salman et al., (2019) Salman, H., Yang, G., Zhang, H., Hsieh, C.-J., and Zhang, P. (2019).  A convex relaxation barrier to tight robustness verification of neural networks.  In Advances in Neural Information Processing Systems, pages 9835–9846.
  * (29) Singh, G., Ganvir, R., Püschel, M., and Vechev, M. (2019a).  Beyond the single neuron convex barrier for neural network certification.  In Advances in Neural Information Processing Systems, pages 15098–15109.
  * Singh et al., (2018) Singh, G., Gehr, T., Mirman, M., Püschel, M., and Vechev, M. (2018).  Fast and effective robustness certification.  In Advances in Neural Information Processing Systems, pages 10802–10813.
  * (31) Singh, G., Gehr, T., Püschel, M., and Vechev, M. (2019b).  An abstract domain for certifying neural networks.  Proceedings of the ACM on Programming Languages, 3(POPL):1–30.
  * Szegedy et al., (2013) Szegedy, C., Zaremba, W., Sutskever, I., Bruna, J., Erhan, D., Goodfellow, I., and Fergus, R. (2013).  Intriguing properties of neural networks.  arXiv preprint arXiv:1312.6199.
  * Tjandraatmadja et al., (2020) Tjandraatmadja, C., Anderson, R., Huchette, J., Ma, W., Patel, K., and Vielma, J. P. (2020).  The convex relaxation barrier, revisited: Tightened single-neuron relaxations for neural network verification.  arXiv preprint arXiv:2006.14076.
  * Tjeng et al., (2018) Tjeng, V., Xiao, K. Y., and Tedrake, R. (2018).  Evaluating robustness of neural networks with mixed integer programming.  In International Conference on Learning Representations.
  * Weng et al., (2018) Weng, L., Zhang, H., Chen, H., Song, Z., Hsieh, C.-J., Daniel, L., Boning, D., and Dhillon, I. (2018).  Towards fast computation of certified robustness for relu networks.  In International Conference on Machine Learning, pages 5276–5285.
  * Wong and Kolter, (2018) Wong, E. and Kolter, Z. (2018).  Provable defenses against adversarial examples via the convex outer adversarial polytope.  In International Conference on Machine Learning, pages 5286–5295. PMLR.
  * Wong et al., (2018) Wong, E., Schmidt, F., Metzen, J. H., and Kolter, J. Z. (2018).  Scaling provable adversarial defenses.  In Advances in Neural Information Processing Systems, pages 8400–8409.
  * Xu et al., (2017) Xu, W., Evans, D., and Qi, Y. (2017).  Feature squeezing: Detecting adversarial examples in deep neural networks.  arXiv preprint arXiv:1704.01155.
  * Zhang et al., (2018) Zhang, H., Weng, T.-W., Chen, P.-Y., Hsieh, C.-J., and Daniel, L. (2018).  Efficient neural network robustness certification with general activation functions.  In Advances in neural information processing systems, pages 4939–4948.

##  5 Appendix

###  5.1 Notation

We introduce the following additional notation:

######  Definition 4.

(Row-wise $q$-norm). Let $\mathbf{A}\in{\mathbb{R}}^{m\times n}$ be an
arbitrary real-valued matrix. Define the row-wise $q$-norm
$||\cdot||_{q\cdot\cdot}$ as the following column vector in
$({\mathbb{R}}^{+}_{0})^{m}$:

$\displaystyle||\mathbf{A}||_{q\cdot\cdot}=(||\mathbf{A}_{1\cdot\cdot}||_{q},||\mathbf{A}_{2\cdot\cdot}||_{q},\dots,||\mathbf{A}_{m\cdot\cdot}||_{q})^{\top}\,.$
(32)

###  5.2 Interval Bound Propagation

Algorithm 1 Interval Bound Propagation

input: Parameters $\\{{\mathbf{W}}^{(i)},{\mathbf{b}}^{(i)}\\}_{i=1}^{L}$,
input ${\mathbf{x}}$, $p$-norm, perturbation size $\epsilon$

${\boldsymbol{\ell}}^{(1)}={\mathbf{W}}^{(1)}{\mathbf{x}}-\epsilon||{\mathbf{W}}^{(1)}||_{p^{*}\cdot\cdot}+{\mathbf{b}}^{(1)}$

${\boldsymbol{u}}^{(1)}={\mathbf{W}}^{(1)}{\mathbf{x}}+\epsilon||{\mathbf{W}}^{(1)}||_{p^{*}\cdot\cdot}+{\mathbf{b}}^{(1)}$

for $i=2,\dots,L$ do

${\boldsymbol{\ell}}^{(i)}={\mathbf{W}}^{(i)}_{+}\sigma({\boldsymbol{\ell}}^{(i-1)})+{\mathbf{W}}^{(i)}_{-}\sigma({\boldsymbol{u}}^{(i-1)})+{\mathbf{b}}^{(i)}$

${\boldsymbol{u}}^{(i)}={\mathbf{W}}^{(i)}_{-}\sigma({\boldsymbol{\ell}}^{(i-1)})+{\mathbf{W}}^{(i)}_{+}\sigma({\boldsymbol{u}}^{(i-1)})+{\mathbf{b}}^{(i)}$

end for

output: bounds
$\\{{\boldsymbol{\ell}}^{(i)},{\boldsymbol{u}}^{(i)}\\}_{i=1}^{L}$

###  5.3 Relaxation based Bound Computation

When each non-linear layer is bounded by exactly one linear lower relaxation
function and one linear upper relaxation function, the relaxed verification
problem can be solved greedily (Salman et al.,, 2019).

In the special case
$\underline{{\boldsymbol{\Sigma}}}^{(i)}=\overline{{\boldsymbol{\Sigma}}}^{(i)}\equiv{\boldsymbol{\Sigma}}^{(i)}$,
we can derive concise closed-form expressions for the lower- and upper pre-
activation bounds. To this end, we replace the activations with their
enveloping relaxation functions
$\underline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})={{\boldsymbol{\Sigma}}}^{(i)}{\mathbf{x}}^{(i)}+\underline{{\boldsymbol{\beta}}}^{(i)}$
resp.
$\overline{{\boldsymbol{\Phi}}}^{(i)}({\mathbf{x}}^{(i)})={{\boldsymbol{\Sigma}}}^{(i)}{\mathbf{x}}^{(i)}+\overline{{\boldsymbol{\beta}}}^{(i)}$
where the offsets
$\underline{{\boldsymbol{\beta}}}^{(i)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(i)}_{j^{\prime}}$
are for each layer $i$ and unit $j^{\prime}$ chosen such that we under-
estimate lower bounds ${\boldsymbol{\ell}}^{(i)}_{j}$ and over-estimate upper
bounds ${\boldsymbol{u}}^{(i)}_{j}$,

$\displaystyle{\boldsymbol{\ell}}^{(i)}_{j}$
$\displaystyle\geq\min_{{\boldsymbol{\xi}}\in\mathcal{B}(\mathbf{0})}\,\min_{{\boldsymbol{\Phi}}^{(i-1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\Phi}}}^{(i-1)}_{j^{\prime}},\overline{{\boldsymbol{\Phi}}}^{(i-1)}_{j^{\prime}}\\}}\cdots\min_{{\boldsymbol{\Phi}}^{(1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\Phi}}}^{(1)}_{j^{\prime}},\overline{{\boldsymbol{\Phi}}}^{(1)}_{j^{\prime}}\\}}\Big{\\{}$
(33)
$\displaystyle{\mathbf{W}}^{(i)}_{j\cdot\cdot}{\boldsymbol{\Phi}}^{(i-1)}(\,\cdots\,{\mathbf{W}}^{(2)}{\boldsymbol{\Phi}}^{(1)}({\mathbf{W}}^{(1)}({\mathbf{x}}+{\boldsymbol{\xi}})+{\mathbf{b}}^{(1)})+{\mathbf{b}}^{(2)}\,\cdots\,)+{\mathbf{b}}^{(L)}\Big{\\}}$
$\displaystyle=\min_{{\boldsymbol{\xi}}\in\mathcal{B}(\mathbf{0})}\,\min_{{\boldsymbol{\beta}}^{(i-1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}}\\}}\cdots\min_{{\boldsymbol{\beta}}^{(1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}}\\}}\Big{\\{}{\mathbf{W}}^{(i)}_{j\cdot\cdot}({\boldsymbol{\Sigma}}^{(i-1)}(\cdots{\mathbf{W}}^{(2)}($
$\displaystyle\quad\quad{\boldsymbol{\Sigma}}^{(1)}({\mathbf{W}}^{(1)}({\mathbf{x}}+{\boldsymbol{\xi}})+{\mathbf{b}}^{(1)})+{{\boldsymbol{\beta}}^{(1)}})+{\mathbf{b}}^{(2)}\cdots)+{{\boldsymbol{\beta}}^{(i-1)}})+{\mathbf{b}}^{(i)}_{j}\Big{\\}}$
$\displaystyle{\boldsymbol{u}}^{(i)}_{j}$
$\displaystyle\leq\max_{{\boldsymbol{\xi}}\in\mathcal{B}(\mathbf{0})}\,\max_{{\boldsymbol{\Phi}}^{(i-1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\Phi}}}^{(i-1)}_{j^{\prime}},\overline{{\boldsymbol{\Phi}}}^{(i-1)}_{j^{\prime}}\\}}\cdots\max_{{\boldsymbol{\Phi}}^{(1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\Phi}}}^{(1)}_{j^{\prime}},\overline{{\boldsymbol{\Phi}}}^{(1)}_{j^{\prime}}\\}}\Big{\\{}$
$\displaystyle{\mathbf{W}}^{(i)}_{j\cdot\cdot}{\boldsymbol{\Phi}}^{(i-1)}(\,\cdots\,{\mathbf{W}}^{(2)}{\boldsymbol{\Phi}}^{(1)}({\mathbf{W}}^{(1)}({\mathbf{x}}+{\boldsymbol{\xi}})+{\mathbf{b}}^{(1)})+{\mathbf{b}}^{(2)}\,\cdots\,)+{\mathbf{b}}^{(L)}\Big{\\}}$
$\displaystyle=\max_{{\boldsymbol{\xi}}\in\mathcal{B}(\mathbf{0})}\,\max_{{\boldsymbol{\beta}}^{(i-1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}}\\}}\cdots\max_{{\boldsymbol{\beta}}^{(1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}}\\}}\Big{\\{}{\mathbf{W}}^{(i)}_{j\cdot\cdot}({\boldsymbol{\Sigma}}^{(i-1)}(\cdots{\mathbf{W}}^{(2)}($
$\displaystyle\quad\quad{\boldsymbol{\Sigma}}^{(1)}({\mathbf{W}}^{(1)}({\mathbf{x}}+{\boldsymbol{\xi}})+{\mathbf{b}}^{(1)})+{{\boldsymbol{\beta}}^{(1)}})+{\mathbf{b}}^{(2)}\cdots)+{{\boldsymbol{\beta}}^{(i-1)}})+{\mathbf{b}}^{(i)}_{j}\Big{\\}}$

Note that the $j^{\prime}$ span different ranges for the different variables.

For the sake of simplicity, we introduce special notation for products over
decreasing index sequences, where the index is counted down from the product-
superscript to the product-subscript

######  Definition 5.

(Products over decreasing index sequences). Let
$\mathbf{A}^{(k)},k\in[m,\dots,n]$ be an arbitrary sequence of (real-valued)
matrices with matching inner dimensions. For ease of notation, define

$\prod^{k=n}_{m}\mathbf{A}^{(k)}=\prod_{k=0}^{n-m}\mathbf{A}^{(n-k)}=\mathbf{A}^{(n)}\mathbf{A}^{(n-1)}\cdots\mathbf{A}^{(m)}\quad\textnormal{where}\,\,\
n\geq m$ (34)

We use the usual convention that the empty product equals one, i.e.
$\prod_{k\in\emptyset}(\cdot)=1$. Thus, $\prod_{k=m}^{n}(\cdot)=1$ and
$\prod^{k=n}_{m}(\cdot)=1$ whenever $n<m$.

We also introduce the following closed-form expressions for the layer-wise
activation-relaxation

######  Definition 6.

(Closed-form expressions for layer-wise activation-relaxation).

$\displaystyle{\boldsymbol{\phi}}^{(i)}$
$\displaystyle=\big{(}\prod^{k=i}_{2}{\mathbf{W}}^{(k)}{\boldsymbol{\Sigma}}^{(k-1)}\big{)}{\mathbf{W}}^{(1)}{\mathbf{x}}+\sum_{j=1}^{i}\big{(}\prod^{k=i}_{j+1}{\mathbf{W}}^{(k)}{\boldsymbol{\Sigma}}^{(k-1)}\big{)}{\mathbf{b}}^{(j)}$
(35) $\displaystyle{\boldsymbol{\Lambda}}^{(j)}$
$\displaystyle=\big{(}\prod^{k=i}_{j+1}{\mathbf{W}}^{(k)}{\boldsymbol{\Sigma}}^{(k-1)}\big{)}{\mathbf{W}}^{(j)}\quad\textnormal{for}\quad
j=1,\dots,i$ (36)

Rewriting the above equations using the explicit expressions for the
concatenation of layers,

$\displaystyle{\boldsymbol{\ell}}^{(i)}_{j}$
$\displaystyle\geq\min_{{\boldsymbol{\xi}}\in\mathcal{B}(\mathbf{0})}\,\min_{{\boldsymbol{\beta}}^{(i-1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}}\\}}\\!\cdots\\!\min_{{\boldsymbol{\beta}}^{(1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}}\\}}\\!\Big{\\{}{\boldsymbol{\phi}}^{(i)}_{j}+{\boldsymbol{\Lambda}}^{(1)}_{j\cdot\cdot}{\boldsymbol{\xi}}+\sum_{k=2}^{i}{\boldsymbol{\Lambda}}^{(k)}_{j\cdot\cdot}{\boldsymbol{\beta}}^{(k-1)}\Big{\\}}$
(37) $\displaystyle{\boldsymbol{u}}^{(i)}_{j}$
$\displaystyle\leq\max_{{\boldsymbol{\xi}}\in\mathcal{B}(\mathbf{0})}\,\max_{{\boldsymbol{\beta}}^{(i-1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(i-1)}_{j^{\prime}}\\}}\\!\cdots\\!\max_{{\boldsymbol{\beta}}^{(1)}_{j^{\prime}}\in\\{\underline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}},\overline{{\boldsymbol{\beta}}}^{(1)}_{j^{\prime}}\\}}\\!\Big{\\{}{\boldsymbol{\phi}}^{(i)}_{j}+{\boldsymbol{\Lambda}}^{(1)}_{j\cdot\cdot}{\boldsymbol{\xi}}+\sum_{k=2}^{i}{\boldsymbol{\Lambda}}^{(k)}_{j\cdot\cdot}{\boldsymbol{\beta}}^{(k-1)}\Big{\\}}$

we can see that the different minimizations and maximizations over optimal
perturbations, lower- and upper- offset vectors decouple, hence we can resolve
them separately, e.g. for $\mathcal{B}_{p}^{\epsilon}(\mathbf{0})$:

$\displaystyle{\boldsymbol{\ell}}^{(i)}$
$\displaystyle\geq{\boldsymbol{\phi}}^{(i)}-\epsilon||{\boldsymbol{\Lambda}}^{(1)}||_{p^{*}\cdot\cdot}+\sum_{j=2}^{i}\big{(}{\boldsymbol{\Lambda}}^{(j)}_{+}\underline{{\boldsymbol{\beta}}}^{(j-1)}+{\boldsymbol{\Lambda}}^{(j)}_{-}\overline{{\boldsymbol{\beta}}}^{(j-1)}\big{)}$
(38) $\displaystyle{\boldsymbol{u}}^{(i)}$
$\displaystyle\leq{\boldsymbol{\phi}}^{(i)}+\epsilon||{\boldsymbol{\Lambda}}^{(1)}||_{p^{*}\cdot\cdot}+\sum_{j=2}^{i}\big{(}{\boldsymbol{\Lambda}}^{(j)}_{-}\underline{{\boldsymbol{\beta}}}^{(j-1)}+{\boldsymbol{\Lambda}}^{(j)}_{+}\overline{{\boldsymbol{\beta}}}^{(j-1)}\big{)}$

where $p^{*}$ is the Hölder conjugate $1/p+1/{p^{*}}=1$ and where
$||\cdot||_{p^{*}\cdot\cdot}$ denotes the row-wise $p^{*}$-norm.

Choosing the same slope for the lower- and upper- relaxation functions
$\alpha^{(i)}_{j}={\boldsymbol{u}}^{(i)}_{j}/({\boldsymbol{u}}^{(i)}_{j}-{\boldsymbol{\ell}}^{(i)}_{j})$,
i.e.
$\underline{{\boldsymbol{\Sigma}}}^{(i)}=\overline{{\boldsymbol{\Sigma}}}^{(i)}$,
recovers Fast-Lin (Weng et al.,, 2018), while adaptively setting
$\alpha^{(i)}_{j}$ to its boundary values depending on which relaxation has
the smaller volume recovers CROWN (Zhang et al.,, 2018).

Algorithm 2 Relaxation based Bound Computation (for
$\underline{{\boldsymbol{\Sigma}}}^{(i)}\equiv\overline{{\boldsymbol{\Sigma}}}^{(i)}$,
$\underline{{\boldsymbol{\beta}}}^{(i)}\equiv\mathbf{0}$)

input: Parameters $\\{{\mathbf{W}}^{(i)},{\mathbf{b}}^{(i)}\\}_{i=1}^{L}$,
input ${\mathbf{x}}$, $p$-norm, perturbation size $\epsilon$

${\boldsymbol{\Lambda}}^{(1)}={\mathbf{W}}^{(1)}$

${\boldsymbol{\phi}}^{(1)}={\mathbf{W}}^{(1)}{\mathbf{x}}+{\mathbf{b}}^{(1)}$

${\boldsymbol{\ell}}^{(1)}={\boldsymbol{\phi}}^{(1)}-\epsilon||{\boldsymbol{\Lambda}}^{(1)}||_{p^{*}\cdot\cdot}$

${\boldsymbol{u}}^{(1)}={\boldsymbol{\phi}}^{(1)}+\epsilon||{\boldsymbol{\Lambda}}^{(1)}||_{p^{*}\cdot\cdot}$

for $i=2,\dots,L$ do

Compute ${\boldsymbol{\Sigma}}^{(i-1)}$,
$\overline{{\boldsymbol{\beta}}}^{(i-1)}$ from ${\boldsymbol{\ell}}^{(i-1)}$,
${\boldsymbol{u}}^{(i-1)}$

${\boldsymbol{\Lambda}}^{(j)}={\mathbf{W}}^{(i)}{\boldsymbol{\Sigma}}^{(i-1)}{\boldsymbol{\Lambda}}^{(j)}$
for $j=1,\dots,i-1$

${\boldsymbol{\Lambda}}^{(i)}={\mathbf{W}}^{(i)}$

${\boldsymbol{\phi}}^{(i)}={\mathbf{W}}^{(i)}{\boldsymbol{\Sigma}}^{(i-1)}{\boldsymbol{\phi}}^{(i-1)}+{\mathbf{b}}^{(i)}$

${\boldsymbol{\ell}}^{(i)}={\boldsymbol{\phi}}^{(i)}-\epsilon||{\boldsymbol{\Lambda}}^{(1)}||_{p^{*}\cdot\cdot}+\sum_{j=2}^{i}\left({\boldsymbol{\Lambda}}^{(j)}_{-}\overline{{\boldsymbol{\beta}}}^{(j-1)}\right)$

${\boldsymbol{u}}^{(i)}={\boldsymbol{\phi}}^{(i)}+\epsilon||{\boldsymbol{\Lambda}}^{(1)}||_{p^{*}\cdot\cdot}+\sum_{j=2}^{i}\left({\boldsymbol{\Lambda}}^{(j)}_{+}\overline{{\boldsymbol{\beta}}}^{(j-1)}\right)$

end for

output: bounds
$\\{{\boldsymbol{\ell}}^{(i)},{\boldsymbol{u}}^{(i)}\\}_{i=1}^{L}$