# Verifying Probabilistic Specifications with Functional Lagrangians

Leonard Berrada  Sumanth Dathathri  Krishnamurthy (Dj) Dvijotham  Robert
Stanforth  Rudy Bunel  Jonathan Uesato  Sven Gowal  M. Pawan Kumar

###### Abstract

We propose a general framework for verifying input-output specifications of
neural networks using _functional Lagrange multipliers_ that generalizes
standard Lagrangian duality. We derive theoretical properties of the
framework, which can handle arbitrary probabilistic specifications, showing
that it provably leads to tight verification when a sufficiently flexible
class of functional multipliers is chosen. With a judicious choice of the
class of functional multipliers, the framework can accommodate desired trade-
offs between tightness and complexity. We demonstrate empirically that the
framework can handle a diverse set of networks, including Bayesian neural
networks with Gaussian posterior approximations, MC-dropout networks, and
verify specifications on adversarial robustness and out-of-distribution (OOD)
detection. Our framework improves upon prior work in some settings and also
generalizes to new stochastic networks and probabilistic specifications, like
distributionally robust OOD detection.

Deep Learning, Probabilistic Inference, Verification

##  1 Introduction

Neural network (NN) verification algorithms seek to prove input-output
specifications, i.e., specifications requiring that the output of the network
satisfies certain constraints for all inputs in a given domain. In this work,
we develop a general framework for verifying specifications of neural networks
composed of an arbitrary sequences of (possibly stochastic) layers. This
allows us to handle both deterministic and stochastic networks (e.g., Bayesian
neural networks (BNNs) or networks using MC-dropout). Furthermore, the
specifications can be defined on output distributions given some input
constraints. This allows us to handle a diverse set of specifications,
including adversarial robustness and uncertainty quantification (See Section
3.1 for concrete examples).

The key to achieving the aforementioned flexibility through our framework is
the generalization of Lagrangian duality, which is often used in verification.
Specifically, instead of using standard linear Lagrange multipliers, we allow
for non-linear functions. We refer to the resulting formulation as _functional
Lagrangian multipliers_. We show that optimizing over sufficiently flexible
classes of functional multipliers is guaranteed to give exact verification
algorithms. We show that even restricted functional multipliers provide a
bound on the specification, which makes the problem computationally tractable.
Our main contributions are:

  * •

We derive a framework for verifying specifications on the probabilistic
predictions of probabilistic neural networks (e.g., BNNs). For example, we can
verify that such a probabilistic neural network never makes high confidence
predictions for a class of OOD inputs.

  * •

We prove that the framework leads to exact verification (i.e. the algorithm
returns verified if and only if the specification is true) if a sufficiently
flexible class of functional multipliers is used. We then show how to easily
achieve trade-offs between computational efficiency and tightness (how far we
are from exact verification) by controlling the multiplier’s complexity.

  * •

We demonstrate empirically that by choosing flexible functional multipliers,
we can improve the tightness of verified bounds on various tasks including
bounding confidence assigned by a classifier given a set of OOD inputs.
Furthermore, we are able to handle new specifications on distributionally
robust OOD detection that previous approaches do not address.

##  2 Related Work

We focus on prior work on verification that addresses probabilistic
specifications. Let $\mathcal{P}_{0}$ be a space of distributions on the
inputs $x$ to a network $\phi$. $\phi\left({p}\right)$ denotes the
distribution over outputs of the network when the inputs follow a distribution
$p$. In this paper, we address the following problem: Verify that: $\forall
p\in\mathcal{P}_{0},\>\phi\left({p}\right)\in\mathcal{P}_{out},$ where
$\mathcal{P}_{out}$ represents a constraint on the output distribution. For
example, a specification can require that the probability of any particular
label is smaller than a threshold when $\mathcal{P}_{0}$ represents a class of
distributions over OOD inputs.

Weng et al. (2019) address certifying robustness of neural networks with
respect to probabilistic perturbations of the input. This corresponds to
choosing a single distribution from $\mathcal{P}_{0}$ rather than finding a
worst case over all distributions in $\mathcal{P}_{0}$. Fazlyab et al. (2019)
address the same problem using semi-definite programming relaxations. However,
these methods do not handle stochastic networks like Bayesian neural networks
or MC dropout networks.

Wicker et al. (2020) defines a notion of probabilistic safety for BNNs:
$P_{w}\left[\forall x\in\mathcal{X},\>\phi_{w}(x)\in\mathcal{C}\right]\geq
p_{\min}$, where $\phi_{w}$ denotes the network with parameters $w$, $P_{w}$
denotes the distribution over network weights (like a Gaussian posterior on
weights, for example), and if the above probability is greater than a desired
threshold ($p_{\min}\in[0,1]$), the network is deemed to be safe. However,
BNNs in practice are used to make probabilistic predictions over the inputs by
averaging over predictions under the distribution of network weights. In this
less restrictive setting, it is sufficient if the constraint is satisfied by
the probabilistic prediction that averages over weights sampled from $P_{w}$:
$\forall x\in\mathcal{X},\>P_{w}\left[\phi_{w}(x)\in\mathcal{C}\right]\geq
p_{\min}$, where $\phi_{w}(x)$ denotes the distribution over outputs for
$x\in\mathcal{X}$. Wicker et al. (2020) observe that
$\min_{x\in\mathcal{X}}P_{w}\left[\phi_{w}(x)\in\mathcal{C}\right]\geq
P_{w}\left[\forall x\in\mathcal{X},\>\phi_{w}(x)\in\mathcal{C}\right]$, making
the second constraint less restrictive.

Our work also generalizes to Bitterwolf et al. (2020b), which studies
specifications of the output distribution of neural networks when the inputs
and network are deterministic. In contrast, our framework’s flexibility allows
for stochastic networks as well. Furthermore, while Bitterwolf et al. (2020b)
are concerned with training networks to be verifiable, their verification
method fails for networks trained in a verification-agnostic manner. In our
experiments, we are able to provide stronger guarantees for networks that are
not explicitly trained to be amenable to verification.

Finally, our framework also generalizes prior work on deterministic neural
network verification based on bound propagation and convex relaxation
(Raghunathan et al., 2018; Zhang et al., 2018; Dvijotham et al., 2018).

##  3 Probabilistic Specifications

##### Problem Formulation:

We consider networks $\phi$ composed of a sequence of possibly stochastic
layers $\pi_{1},\ldots\pi_{K}$, where $\pi_{k}\left({y|x}\right)$ denotes the
probability that the output of layer $k$ is $y$ when its input is $x$. We
write $y\sim\pi_{k}\left({x}\right)$ to denote that $y$ is sampled from the
distribution over outputs of layer $k$ given input $x$ to layer $k$.

Let $p_{0}$ denote a distribution over the inputs $\mathcal{X}$. We would like
to verify some property of the distribution of outputs of the model when the
input is sampled from $p_{0}$ (for example, a Gaussian centered at some
nominal input $x_{0}$). Unless otherwise specified, we assume that the set of
outputs $\mathcal{Y}$ is a subset of $\mathbb{R}^{l}$, where $l$ is the number
of labels, and the output of the network are logits corresponding to
unnormalized log-confidence scores assigned to the labels $\\{1,\ldots,l\\}$.

######  Definition 1 (Probabilistic verification problem).

Given a (possibly stochastic) neural network
$\phi:\mathcal{X}\to\mathcal{P}\left({\mathcal{Y}}\right)$, a set of
distributions over the input $\mathcal{P}_{0}$ and a functional
$\boldsymbol{\psi}:\mathcal{P}\left({\mathcal{Y}}\right)\mapsto\mathbb{R}$,
the probabilistic verification problem is to check that the following is true:

$\displaystyle\forall\>p_{0}\in\mathcal{P}_{0},\>\boldsymbol{\psi}\left({\phi\left({p_{0}}\right)}\right)\leq
0.$ (1)

We consider the following optimization version:

$\displaystyle\texttt{OPT}=\max_{p_{0}\in\mathcal{P}_{0}}\boldsymbol{\psi}\left({\phi\left({p_{0}}\right)}\right),$
(2)

Having $\texttt{OPT}\leq 0$ here is equivalent to satisfying specification (1)
.

##### Stochastic Networks:

We consider stochastic networks with layers of the form
$\pi_{k}\left({x}\right)=\mathrm{\texttt{ReLU}}\left({\tilde{w}x+\tilde{b}}\right)$,
where $\tilde{w}$ and $\tilde{b}$ are random variables (following a diagonal
Gaussian distribution for BNNs and a Bernoulli-like distribution for MC-
Dropout) - both Bayesian Neural Networks (BNNs) and MC-Dropout Networks are
special cases of this.

We define all specifications on the expected softmax output of the network. In
practice, inference in stochastic networks is done with a finite set of
independent and identically distributed (i.i.d) realizations of the
stochasticity in the layers, which is related to the expected softmax value by
standard concentration of measure arguments. This stochasticity in a layer
$\pi_{k}$ is assumed to be statistically independent of the stochasticity at
other layers.

###  3.1 Examples of Specifications

##### Robustness to Perturbations:

Our goal is to evaluate whether a neural network is robust to norm-bounded
adversarial perturbations. To model this within our framework, we define as a
set of $\delta$ input distributions centered at points within an $\epsilon$
ball of a nominal point $x_{\text{nom}}$, with label $i\in\\{1,\ldots,l\\}$:
$\mathcal{P}_{0}=\\{\delta_{x}:\left\|{x-x_{\text{nom}}}\right\|\leq\epsilon\\}$.
We define the specification to be
$\boldsymbol{\psi}\left({p}\right)=\operatorname*{\mathbb{E}}_{y\sim
p}\left[{y_{i}-y_{j}}\right]$, the expected difference between the logits
assigned to the labels $j$ and $i$, and an adversary fails if the optimal
value of (2) is smaller than $0$ for this $\boldsymbol{\psi}\left({.}\right)$.
If the network and the output are deterministic, (2) corresponds to optimizing
the difference between the logits for labels $j$ and $i$.

##### OOD Detection for Stochastic Neural Networks:

Here, we consider the problem of verifying that a stochastic neural network
assign low confidence scores to all labels for inputs that are OOD. We
consider two formulations here:

_Robust OOD detection:_ Bitterwolf et al. (2020a) consider classes of inputs
that are out of distribution. For example, those generated by perturbations of
a nominal OOD point:
$\mathcal{P}_{0}=\\{\delta_{x}:\left\|{x-x_{ood}}\right\|\leq\epsilon\\}$. We
would like to verify that the maximum confidence assigned by the classifier to
any of these points is small. To this end, we define:

$\displaystyle\boldsymbol{\psi}\left({p}\right)=\operatorname*{\mathbb{E}}_{y\sim
p}\left[{\frac{\exp\left({y_{i}}\right)}{\sum_{j}\exp\left({y_{j}}\right)}}\right],$
(3)

and we require that (2) is smaller than a maximum confidence threshold
$p_{\max}$ for each $i$.

_Distributionally robust OOD detection:_ The above specification assumes that
$\mathcal{P}_{0}$ is a set of $\delta$ distributions. However, there are often
situations where we would like to have a whole class of OOD distributions be
classified (on average) with low confidence. This can be expressed as
$\mathcal{P}_{0}=\\{x_{ood}+\omega:\omega\in\mathcal{P}_{noise}\\}$, where
$\mathcal{P}_{noise}$ is a class of noise distributions. We define the
expected softmax as in (3) and require that (2) is smaller than a maximum
confidence threshold $p_{\max}$ for each $i$.

##  4 Functional Lagrangians for Verification

Problem (2) is a challenging nonlinear and possibly stochastic optimization
problem. Thus, solving it to global optimality in general is computationally
intractable. However, in order to verify the specification, it suffices to
compute an upper bound on the optimal value. Lagrangian relaxation is a
standard technique from constrained optimization to do this in the
deterministic setting and it has indeed been successfully used in neural
network verification (Dvijotham et al., 2018). For the sake of clarity, we
start by introducing this technique and motivate our framework.

##### Problem Setup:

Consider a deterministic neural network composed of a series of layers
$h_{1},\ldots,h_{K}$. Denote the activations received as output of layer
$h_{k}$ as $x_{k}\in\mathbb{R}^{n_{k}}$ and the input to the network as
$x_{0}\in\mathbb{R}^{n_{0}}$ .

In this setting, problem (2) can be written as:

$\texttt{OPT}=\max_{x_{0}\in\mathcal{X}_{0}}\psi\circ h_{K}\circ\ldots
h_{1}\left({x_{0}}\right),$

where $\circ$ denotes function composition. This problem can be rephrased as a
constrained optimization problem as follows:

$\displaystyle\texttt{OPT}=$
$\displaystyle\max_{x_{0},\ldots,x_{K}}\psi\left({x_{K}}\right)$ (4) s.t.:
$\displaystyle
x_{0}\in\mathcal{X}_{0},\>x_{k}=h_{k}\left({x_{k-1}}\right)\forall\>k\in\\{1,\ldots,K\\},$

where $\mathcal{X}_{0}$ denotes the set of inputs over which the constraint
must be satisfied (for example, an $\epsilon$-ball around an input).

##### Lagrangian Relaxation:

The traditional Lagrangian relaxation (Bertsekas, 2015) of this problem is:

$\max_{\mathbf{x}}\psi\left({x_{K}}\right)+\sum_{k=1}^{K}\theta_{k}^{T}\left({h_{k}\left({x_{k-1}}\right)-x_{k}}\right),$
(5)

where the Lagrange multipliers $\theta_{k}\in\mathbb{R}^{k}$ can be chosen
arbitrarily and $\mathbf{x}=\\{x_{k}\\}_{k=0}^{K}$ with
$x_{0}\in\mathcal{X}_{0}$. Thus, we have converted a constrained optimization
problem to a simpler unconstrained optimization problem (or problem with only
simple domain constraints).

Consider any $\mathbf{x}$ that satisfies the constraints of (4). Then, the
objective functions of the two problems coincide in this case since all the
terms involving $\theta_{k}$ are multiplied by a $0$ vector (since we have
$x_{k}=h_{k}\left({x_{k-1}}\right)$). Thus, the optimal value of (5) is an
upper bound on OPT for any choice of
$\boldsymbol{\theta}=\\{\theta_{k}\\}_{k=0}^{K-1}$. This property is known as
weak duality (Bertsekas, 2015) . Thus, if we can find a choice of
$\boldsymbol{\theta}$ such that the optimal objective value of (5) is smaller
than $0$, then necessarily $\texttt{OPT}\leq 0$ and the specification is
verified. However, since the upper bound is not tight, the converse is not
true. In other words, even if there is no such value of $\boldsymbol{\theta}$,
it is possible that the specification is still true.

##### Limitations of Lagrangian Relaxation:

Given the stated goal of solving problem (2) as accurately as possible, there
are two limitations in the Lagrangian relaxation approach. First, as explained
above, since the relaxation is not tight in general, the resulting upper bound
can be too conservative in practice. It would therefore be desirable to have a
more flexible approach that can tighten the relaxation in order to solve
difficult problems more accurately. Second, the approach so far does not take
into account the stochasticity of the problem as formulated in (2).

##### Functional Lagrangians:

In what follows, we describe how our framework alleviates these two issues. We
describe the functional Lagrangians for the deterministic setting here to
convey the key intuition and state the general framework that applies to the
stochastic setting in section 4.1. The objective in (5) is a linear function
of the violation in the constraints of (4), i.e.,
$x_{k}-h_{k}\left({x_{k-1}}\right)$. However, it is easy to see that weak
duality holds even if we modify (5) to have a nonlinear term of the form
$\lambda_{k}\left({x_{k}}\right)-\lambda_{k}\left({h_{k}\left({x_{k-1}}\right)}\right)$:

$\begin{split}&\max_{x_{0}\in\mathcal{X}_{0},\ldots,x_{K}}\psi\left({x_{K}}\right)+\sum_{k=1}^{K}\lambda_{k}\left({h_{k}\left({x_{k-1}}\right)}\right)-\lambda_{k}\left({x_{k}}\right)\\\
&\quad=\max_{x_{K}}\boldsymbol{\psi}\left({x_{K}}\right)-\lambda_{K}\left({x_{K}}\right)\\\
&\qquad+\sum_{k=1}^{K-1}\max_{x_{k}}\lambda_{k+1}\left({h_{k+1}\left({x_{k}}\right)}\right)-\lambda_{k}\left({x_{k}}\right)\\\
&\qquad+\max_{x_{0}\in\mathcal{X}_{0}}\lambda_{1}\left({h_{0}\left({x_{0}}\right)}\right),\end{split}$
(6)

where $\lambda_{k}:\mathbb{R}^{n_{k}}\to\mathbb{R}$ is an arbitrary real-
valued function. If we consider $\lambda_{k}\left({x}\right)=\theta_{k}^{T}x$,
we recover the traditional Lagrangian. However, by allowing $\lambda_{k}$ to
be nonlinear, we can derive better bounds on the objective (5).

###  4.1 The Functional Lagrangian Framework

Having described the concept of functional Langrangians, we are now ready to
deploy it for verification.

##### Assumptions:

In order to derive our verification framework, we make the following
assumptions on the activations at each layer:
_Assumption 1_ : $\exists l_{0}\leq u_{0}\in\mathbb{R}^{n}$ such that for each
feasible input distribution $p_{0}\in\mathcal{P}_{0}$,
$\text{Support}\left({p_{0}}\right)\subseteq\mathcal{X}_{0}=[l_{0},u_{0}]$.
_Assumption 2_ : Each layer is such that if
$x\in\mathcal{X}_{k}=[l_{k},u_{k}]$, then
$\text{Support}\left({\pi_{k}\left({x}\right)}\right)\subseteq\mathcal{X}_{k+1}=[l_{k+1},u_{k+1}]$.

Assumption 1 is natural since the inputs to neural networks are bounded.
Assumption 2 is somewhat restrictive: it requires that the layer output is
bounded with probability 1, which is not true, for example, if we have a
Bayesian neural network with a Gaussian posterior. However, we can relax this
assumption to requiring that the output is bounded with high probability, as
in Wicker et al. (2020).

Let $\mathcal{X}_{k}$ denote the space of activations at layer $k$. Let
$p_{k}$ denote the distribution of activations at layer $k$ when the inputs
follow distribution $p_{0}$ (so that $p_{K}=\phi\left({p_{0}}\right)$).

##### Functional Lagrangian Dual:

We define the functional Lagrangian dual function as:

$\displaystyle
g\left({\lambda}\right)=\max_{p_{0}\in\mathcal{P}_{0}}g_{0}(p_{0},\lambda_{1})+\sum_{k=1}^{K-1}\max_{x_{k}\in\mathcal{X}_{k}}g_{k}(x_{k},\lambda_{k},\lambda_{k+1})$
$\displaystyle\qquad\quad+\max_{p_{K}\in\mathcal{P}_{K}}g_{K}(p_{K},\lambda_{K}),$
(7) $\displaystyle\text{where:
}g_{0}(p_{0},\lambda_{1})=\operatorname*{\mathbb{E}}_{{x\sim
p_{0}}}\left[{\operatorname*{\mathbb{E}}_{y\sim\pi_{0}\left({x}\right)}\left[{\lambda_{1}\left({y}\right)}\right]}\right],$
$\displaystyle
g_{k}(x_{k},\lambda_{k},\lambda_{k+1})=\operatorname*{\mathbb{E}}_{y\sim\pi_{k}\left({x_{k}}\right)}\left[{\lambda_{k+1}\left({y}\right)}\right]-\lambda_{k}\left({x_{k}}\right),$
$\displaystyle
g_{K}(p_{K},\lambda_{K})=\boldsymbol{\psi}\left({p_{K}}\right)-\operatorname*{\mathbb{E}}_{y\sim
p_{K}}\left[{\lambda_{K}\left({y}\right)}\right].$

The dual $g\left({\lambda}\right)$ can be seen as a generalization of the
Lagrangian relaxation (5) with the following two main modifications: (i) layer
outputs are now described by expectations, and (ii) Lagrangian penalties are
expressed as arbitrary functions $\lambda_{k}\left({x}\right)$ instead of
being restricted to linear functions.

##### Main Result:

Our main result relates the functional Lagrangian dual to the specification
(2).

######  Theorem 1.

For any collection of functions
$\lambda=(\lambda_{1},\ldots,\lambda_{K})\in\mathbb{R}^{\mathcal{X}_{1}}\times\ldots\times\mathbb{R}^{\mathcal{X}_{K}}$,
we have that $g\left({\lambda}\right)\geq\texttt{OPT}$. In particular, if a
choice of $\lambda$ can be found such that $g\left({\lambda}\right)\leq 0$,
then specification (1) is true. Further, if
$\boldsymbol{\psi}\left({p_{K}}\right)=\operatorname*{\mathbb{E}}_{y\sim
p_{K}}\left[{c\left({y}\right)}\right]$, and we define $\lambda^{\star}$ as:

$\displaystyle\lambda^{\star}_{K}\left({x}\right)$
$\displaystyle=c\left({x}\right),$
$\displaystyle\forall\>k\in\\{K-1,\ldots,1\\},\>\lambda^{\star}_{k}\left({x}\right)$
$\displaystyle=\operatorname*{\mathbb{E}}_{y\sim\pi_{k}\left({x}\right)}\left[{\lambda^{\star}_{k+1}\left({y}\right)}\right],$

then $g\left({\lambda^{\star}}\right)=\texttt{OPT}$.

###### Proof.

We give a brief sketch of the proof - the details are in Appendix A. The
layers of the network can be expressed as linear constraints between the
probability distributions over activations at layer $k$ and $k+1$:

$p_{k+1}\left({y}\right)=\int
p_{k}\left({x}\right)\pi_{k}\left({y|x}\right)dx.$

Thus, (2) can be rewritten as an infinite dimensional optimization with
decision variables $p_{0},p_{1},\ldots,p_{K}$ and linear constraints relating
$p_{k}$ and $p_{k+1}$. The Lagrangian dual of this optimization problem has
objective $g\left({\lambda}\right)$ where $\lambda\left({y}\right)$ represent
the Lagrange multipliers for the linear constraint written above (we have
independent $\lambda\left({y}\right)$ for each $y$, hence $\lambda$ is a
_functional multiplier_). By weak duality, we have the result. The second part
of the theorem is easily observed by plugging in $\lambda^{\star}$ in
$g\left({\lambda}\right)$ and observing that the resulting optimization
problem is equivalent to (2). ∎

##### Relationship to Prior Work:

In Appendix F, we show that for feedforward networks, the Functional
Lagrangian Dual with linear multipliers is equivalent to the Lagrangian Dual
from Dvijotham et al. (2018). This is also the dual of the LP relaxation
(Ehlers, 2017) and the basis for most of the efficient neural network
verification algorithms used till date ((Zhang et al., 2018; Singh et al.,
2018), for example) as shown in Liu et al. (2019).

We also prove that with quadratic multipliers and a particular grouping of
layers, our framework is equivalent to the Lagrangian dual of the SDP
formulation from Raghunathan et al. (2018).

Finally, we note that similar mathematical ideas on nonlinear Lagrangians have
been explored in the optimization literature (Nedich & Ozdaglar, 2008;
Feizollahi et al., 2017) but only as a theoretical construct - this literature
has not lead to practical algorithms and do not exploit the staged structure
of optimization problems arising in neural network verification. Further,
these approaches do not handle stochastic networks.

##  5 Practical Verification Algorithms

The choice of functional multipliers affects the difficulty of evaluating
$g\left({\lambda}\right)$, and hence we cannot arbitrarily choose functional
Lagrangians. In fact, since neural network verification is NP-hard (Katz et
al., 2017), we know that computing $g\left({\lambda^{\star}}\right)$ is
intractable in the general case. We need to restrict the form of the
functional Lagrangian to allow for efficient computation of the optimal value
(or upper bound) of $g\left({\lambda}\right)$. We describe several such
choices below:

##### Parameterization:

We now present some practical verification algorithms that instantiate the
functional Lagrangian framework for specific parameterized classes of
Lagrangian functions $\lambda$. We denote the parameterized $\lambda$
functions as
$\lambda\left({\theta}\right)=\left\\{\lambda_{k}\left({x}\right)=\lambda_{k}\left({x;\theta_{k}}\right)\right\\}_{k=1}^{K}$.
Choosing the right class of functions $\lambda(\theta)$ is a trade-off: for
very simple classes (such as linear functions), $g(\lambda(\theta))$ is easy
to compute but may be a loose upper bound on (2), while more expressive
choices lead to tighter relaxation of (2) at the cost of more difficult
evaluation (or bounding) of $g(\lambda(\theta))$.

With some abuse of notation, for convenience, we write
$g_{0}\left({x_{0},\lambda_{0},\lambda_{1}}\right)\coloneqq
g_{0}\left({p_{0},\lambda_{1}}\right)$ and
$g_{K}\left({x_{K},\lambda_{K},\lambda_{K+1}}\right)\coloneqq
g_{K}\left({p_{K},\lambda_{K}}\right)$, with $\lambda_{0}=\lambda_{K+1}=0$.
Then the problem of obtaining the best bound can be written as:

$\min_{\theta}\sum_{k=0}^{K}\max_{x_{k}}g_{k}\left({x_{k},\lambda_{k},\lambda_{k+1}}\right),$
(8)

where the inner maximizations are understood to be performed over the
appropriate domains ($\mathcal{P}_{0}$ for $x_{0}$, $\mathcal{X}_{k}$ for
$x_{k}$, $l=1,\ldots,K-1$ and $\mathcal{P}_{K}$ for $x_{K}$). The overall
procedure is described in Algorithm 1: $\theta$ is minimized by a gradient-
based method in the outer loop; in the inner loop, the decomposed maximization
problems over the $x_{k}$ get solved, potentially in parallel. These inner
problems need not to be solved exactly as long as they can provide sufficient
information about the outer descent direction for $\theta$. We consider three
options depending on the choice of multiplier and type of later (details in
Section 5.1):
1) Solve the maximum in closed-form.
2) Compute a closed-form upper bound on the maximum.
3) Compute an approximate maximum using a local search procedure like gradient
ascent.

To state the algorithm succinctly, we denote this (possibly approximate)
maximization as
$\displaystyle{\max_{x_{k}}}^{\ast}g_{k}\left({x_{k},\lambda_{k},\lambda_{k+1}}\right)$.

For the final verification certificate to be valid, we do require the final
evaluation to provide the exact value of $g(\lambda(\theta^{(T)}))$ (or an
upper bound). In addition, if at any point during the minimization over
$\theta$ we know that $g\left({\lambda\left({\theta}\right)}\right)\leq 0$, we
can stop the algorithm since then (1) is guaranteed to be true.

Algorithm 1 Verification with Functional Lagrangians

Input: initial dual parameters $\theta^{(0)}$, learning-rate $\eta$, number of
iterations $T$.

// Functional Multiplier Training

for $t=0,\ldots,T-1$ do

// potentially in parallel:

for $k=0$ to $K$ do

$d_{\theta}^{(k)}=\nabla_{\theta}\displaystyle{\max_{x_{k}}}^{\ast}g_{k}(x_{k},\lambda_{k},\lambda_{k+1})$

end for

$d_{\theta}=\sum_{k=0}^{K}d_{\theta}^{(k)}$

$\theta^{(t+1)}=\theta^{(t)}-\eta d_{\theta}$.

// or any gradient based optimization

end for

//Functional Multiplier Eval

Return: Exact evaluation of $g(\lambda(\theta^{(T)}))$ (or upper bound)

###  5.1 Solving Inner Maximization Problems

The following paragraphs describe how we solve inner maximization problems.
For neural network (NN) multipliers and linear/quadratic multipliers, we use
methods that were developed in prior work - these are summarized in Table 1.
For the input and output layer and exponential multipliers, we need to develop
custom methods that we describe here.

Layer type | Multiplier | Inner Max (Train) | Inner Max (Eval)
---|---|---|---
Deterministic | NN | PGA | MIP
Stochastic | Quadratic | FO-SDP | FO-SDP
Stochastic | Linear | Closed-Form | Closed-Form

Table 1: Methods for solving inner maximization problems at train/eval time.
MIP=Mixed Integer Program formulation (Tjeng et al., 2017), FO-SDP=First-order
Semidefinite Programming (Dathathri et al., 2020), PGA=Projected Gradient
Ascent

##### NN Multipliers:

Suppose that $\lambda_{k}$ are themselves neural networks with a single real-
valued output and that the network being verified is deterministic. Then we
solve the inner maximization problems approximately via projected gradient
methods at training time, and exactly with a mixed integer program (MIP)
(Tjeng et al., 2017) at evaluation time. Since the MIP formulation is only
valid for ReLU networks, we use ReLU activations. We find that regularizing
the multipliers $\lambda$ during training speeds up the MIP solution during
evaluation. It also helps projected gradient methods solve the maximization
more reliably (See Appendix B.4 for details).

##### Quadratic Multipliers:

If the networks are stochastic, we need to be able to compute the expectations
$\operatorname*{\mathbb{E}}_{z\sim\pi_{k}}\left[{\lambda_{k+1}\left({z}\right)}\right]$.
This can be done easily if the layer is of the form $ws\left({x}\right)+b$,
where $s$ is an element-wise activation function and $(w,b)$ follow a
distribution whose first and second moments can be easily computed. The
resulting inner maximization over $g_{k}$ can be relaxed to an SDP. Following
the approach of Dathathri et al. (2020), the SDP can be rewritten as an
eigenvalue optimization problem with additional inner dual variables
corresponding to the relaxation of the inner maximization. We append these
inner dual variables to $\theta$ and jointly optimize the parameters of the
functional Lagrangian and these additional dual variables in Algorithm 1. In
this case, no separate evaluation phase is needed since the SDP relaxation
provides an upper bound on $\max g_{k}$. Details are in Appendix B.5. For the
special case of linear multipliers, we can bypass the SDP and solve the inner
maximization in closed form, both for deterministic and stochastic networks
(See Appendix B.6).

##### Input Layer Maximization with Linexp Multipliers:

We consider the setting where we verify a deterministic network with
stochastic inputs and constraints on the input distribution
$p_{0}\in\mathcal{P}_{0}$. In this setting, the functional Lagrangian dual
$g\left({\lambda}\right)$ only depends on the input distribution via $g_{0}$,
and $g_{0}$ evaluates to:

$g_{0}\left({p_{0},\lambda_{1}}\right)=\operatorname*{\mathbb{E}}_{x\sim
p_{0}}\left[{\lambda_{1}\left({x}\right)}\right].$

If we choose $\lambda_{1}$ to be a linear or quadratic function, then
$g\left({\lambda}\right)$ only depends on the first and second moments of
$p_{0}$. This implies that the verification results will be unnecessarily
conservative as they don’t use the full information about the distribution
$p_{0}$. To consider the full distribution it suffices to add an exponential
term which evaluates to the moment generating function of the input
distribution. In particular, we consider:

$\mathcal{P}_{0}=\\{\mu+\omega:\omega\sim\mathcal{P}_{noise}\\},$

where $\mathcal{P}_{noise}$ denotes a class of $0$ mean noise distributions
that all satisfy the property of having sub-Gaussian tails (this is true for
many common noise distributions including Bernoulli, Gaussian, truncated
Gaussian):

$\operatorname*{\mathbb{E}}\left[{\exp\left({t\omega_{i}}\right)}\right]\leq\exp\left({\frac{t^{2}\sigma^{2}}{2}}\right)\quad\forall
i,\forall t\in\mathbb{R}.$

We also assume that each component of the noise $\omega_{i}$ is i.i.d. Let
$\lambda_{1}\left({x}\right)=\alpha^{T}x+\exp\left({\gamma^{T}x+\kappa}\right)$,
$\lambda_{2}\left({x}\right)=\beta^{T}x$,
$g_{0}^{\star}=\max_{p_{0}\in\mathcal{P}_{0}}g_{0}\left({p_{0},\lambda_{1}}\right)$,
and
$g_{1}^{\star}=\max_{x\in\mathcal{X}_{1}}g_{0}\left({x,\lambda_{1},\lambda_{2}}\right)$.

######  Proposition 1.

In the setting described above, and with $s$ as the element-wise activation
function:

$\displaystyle
g_{0}^{\star}\leq\alpha^{T}\left({w\mu+b}\right)+\exp\left({\frac{\left\|{w^{T}\gamma}\right\|^{2}\sigma^{2}}{2}+\gamma^{T}b+\kappa}\right),$
$\displaystyle
g_{1}^{\star}\leq\max_{\begin{subarray}{c}x\in\mathcal{X}_{2}\\\
z=s\left({x}\right)\end{subarray}}\beta^{T}\left({w_{2}z+b_{2}}\right)-\alpha^{T}x-\exp\left({\gamma^{T}x+\kappa}\right).$

The maximization in the second equation can be bounded by solving a convex
optimization problem (Appendix B.3).

##### Output Layer Problem:

We focus on specifications that are expressed in terms of an expected value:
$\boldsymbol{\psi}\left({p}\right)=\operatorname*{\mathbb{E}}_{x\sim
p}\left[{c\left({x}\right)}\right]$. In this case, the final layer
maximization can be expressed as:

$\displaystyle\max_{x\in\mathcal{X}_{K}}c\left({x}\right)-\lambda_{K}\left({x}\right).$
(9)

We consider $\lambda_{K}$ that are diagonal quadratic functions:
$\lambda_{K}\left({x}\right)=\alpha^{T}x+\beta^{T}\left({x\odot x}\right)$.
and allow $c\left({x}\right)$ to be an arbitrary linear function of the
softmax function, so we can enforce specifications that are expressible as a
linear function of the softmax probabilities produced by the network:
$c\left({x}\right)=\mu^{T}\text{softmax}\left({x}\right)=\frac{\sum_{i}\mu_{i}\exp\left({x_{i}}\right)}{\sum_{i}\exp\left({x_{i}}\right)}$.

However, problem (9) is not concave, and it is necessary to find the exact
maximum of problem (9) or a guaranteed upper bound for valid verification. The
following result alleviates this by providing an efficient algorithm to bound
(9) for affine $\lambda_{K}$ (more general results in Appendix B.2).

######  Proposition 2.

Assume that $\lambda_{K}$ is an affine function and that $\mu$ is a one-hot
encoding vector. Let $[t_{1},t_{2},\ldots,t_{N}]$ be an ordered sequence of
real numbers such that $t_{1}=\min_{x\in\mathcal{X}_{k}}c(x)$ and
$t_{N}=\max_{x\in\mathcal{X}_{k}}c(x)$. For $i=1,\ldots,N-1$, let $M_{i}$ be:

$\displaystyle\max_{x\in\mathcal{X}_{K}}-\lambda_{K}(x)\text{ subject to:
}t_{i}(1^{\top}\exp(x-\mu^{\top}x))\leq 1$

Then $\max_{1\leq i\leq
N-1}M_{i}+t_{i+1}\geq\max_{x\in\mathcal{X}_{K}}c\left({x}\right)-\lambda_{K}\left({x}\right)$.

Each evaluation of $M_{i}$ can be done efficiently since the objective is
linear and the constraints are convex.

Additional methods to solve (9) are provided in appendix B.2. Namely we show
in section B.2.2 how to exactly solve problem (9) for affine multipliers, at
the cost of a time complexity scaling exponentially with the problem
dimension. In practice, the algorithm runs fast in our experiments since the
dimension (number of classes) is typically only 10. In addition, section B.2.3
describes a tractable upper bound on problem (9) when $\lambda_{K}$ is
quadratic (instead of affine as assumed in Proposition 2).

##  6 Experiments

This sections aims at empirically validating the theoretical flexibility of
the framework and its applicability to both deterministic and stochastic
networks across several specifications. In particular, we show that (1) by
increasing the complexity of the functional multipliers, one can obtain
increasingly tighter bounds, and (2) the functional Lagrangian framework
permits verification of probabilistic neural networks by effectively handling
parameter stochasticity.

The empirical findings are organized in four subsections. First, we verify
standard deterministic convolution neural networks (CNNs) and multi-layer
perceptrons (MLP) for adversarial classification robustness. Second, we tackle
robust OOD detection on deterministic networks. Third, we consider Bayesian
neural networks on similar OOD robustness tasks. Lastly, we tackle robust OOD
detection on for neural networks with noisy perturbations of the input.

###  6.1 Adversarial Robustness on Deterministic NNs

##### Verification Task:

We consider PGD-NN – a 3-layer MLP trained using adversarial training (Madry
et al., 2017) on MNIST (LeCun & Cortes, 2010) – from Raghunathan et al. (2018)
and verify robustness to perturbations of $\ell_{\infty}$ magnitude
$\epsilon=.1$ (when the images are rescaled to $[0,1]$).

Multiplier |  | MLP | MLP |  |  |
---|---|---|---|---|---|---
| Linear | (5 units) | (20 units) | Quad | Best | SDP-FO
Runtime (in $10^{3}$ s) | 25.0 | 43.1 | 123.0 | 43.0 | – |  $\sim$1
Upper Bound | -3.71 | -4.05 | -4.35 | -4.55 | -4.73 | -6.14
Verified Robust Acc (%) | 72 | 76 | 77 | 79 | 79 | 88

Table 2: Runtime, upper bound (averaged across all incorrect labels) and the
verified robust accuracy with various multipliers, and SDP-FO (Dathathri et
al., 2020). The results reported are the average across first 100 samples from
the test subset of the MNIST dataset. As the multipliers grow in complexity,
the runtime increases while the relaxation becomes tighter. _Best_ corresponds
to the best bound obtained across all choices of multipliers. Runtimes are
measured on a machine with a P100 GPU, 32GB RAM.

Training | Perturbation | OOD | GAUC (Lower Bound) | AAUC
---|---|---|---|---
Data | Radius | Dataset | LP-BP | Ours | Upper Bound
MNIST | $0.01$ | EMNIST | 0.65 | 0.79 | 0.91
$0.03$ | (Letters) | 0.10 | 0.15 | 0.86
SVHN | $0.01$ | CIFAR-10 | 0.65 | 0.86 | 0.95
$0.03$ | 0.01 | 0.08 | 0.81
CIFAR-10 | $\frac{1}{255}$ | SVHN | 0.30 | 0.40 | 0.75

Table 3: Area Under Curver (AUC) values: comparing guaranteed AUC (GAUC)
values (verified lower bound on the true AUC values) from LP-BP (baseline) and
our approach with the adversarial AUC (upper bound on the true AUC). Our
approach results in higher guaranteed AUC scores across settings.

##### Method:

We compute bounds on all the intermediate activations using the approach from
Bunel et al. (2020), which proposes an efficient solver for the LP relaxations
of neural networks developed in Ehlers (2017). Using these intermediate
bounds, for every feasible incorrect target label, we compute an upper bound
on the difference between the logits corresponding to the incorrect label and
the true label using (8), and we report that a sample is verified to be
adversarially robust if a non-positive upper bound can be verified for each
target label.

Starting with the LP relaxation (Ehlers, 2017) (which is equivalent to using
our approach with linear multipliers, see Proposition 7), our primary
objective here is to demonstrate that as we increase the complexity of the
multipliers (to non-linear MLPs and quadratic multipliers), we obtain
progressively tighter relaxations. For MLP multipliers, we use the
optimization approach described in Appendix B.4, and for quadratic
multipliers, see Appendix B.5.

##### Results:

The results are summarized in Table 2 (averaged across first 100 samples from
the test subset of the MNIST dataset). We obtain systematically improving
bounds as the complexity of the multiplier increases from linear to MLPs with
an increasing number of ReLU units. In Table 2, since the SCIP solver (Gamrath
et al., 2020a, b) does not always terminate, we use a cascaded approach to
computing the bounds and accuracy, where for a given multiplier we use the
best upper bound from across all multipliers of lesser complexity (See
Appendix C for details). In Figure 2 (Appendix C), we compare the upper bounds
computer from the different multipliers (without cascading) for 50 problem
instances where all approaches terminate. In Table 2, we compare the runtimes
for verification with the various multipliers. We see that while we obtain
tighter relaxations with growing complexity of the multipliers, we also see an
increase in the computation cost marking a trade off between the computation
cost and the tightness of the relaxation. While these results fall short of
SOTA results (SDP-FO) - we believe that further work on scaling the
optimization of $g\left({\lambda}\right)$ is necessary to outperform the SDP-
FO method and is an exciting direction for future work. While the upper bound
improves for ‘Best’ when compared to those obtained with MLPs/quadratic
multipliers, the accuracy does not improve as the improvement in bounds
correspond to instances that are already verified/are far from being verified.

###  6.2 Robust OOD Detection on Deterministic NNs

##### Verification Task:

We bound the largest softmax probability across all labels for out-of-
distribution inputs over $\ell_{\infty}$ perturbations of the input. We train
a model using the CEDA method (Bitterwolf et al., 2020a) with i) MNIST
(handwritten digits) ii) CIFAR-10 (Krizhevsky, 2009) and iii) SVHN (Netzer et
al., 2011) as the in-distribution datasets, and i) a subset of EMNIST (Cohen
et al., 2017) corresponding to English letters, ii) SVHN, and iii) CIFAR as
the corresponding OOD datasets respectively. On the first 1000 samples of the
test subset of the OOD datasets, for the trained networks we seek to verify
that they does not assign high confidence predictions to OOD inputs.

##### Method:

In Bitterwolf et al. (2020a), the authors develop a verification method based
on bound propagation and use this to train networks such that they perform
verifiably robust OOD detection. In contrast, we consider the more challenging
setting of verifying the ability of networks to perform OOD detection, where
the networks are trained in a manner agnostic to verification. For both, the
bound propagation based baseline (LP-BP) from (Bitterwolf et al., 2020a) and
our approach (with linear functional multipliers), we use bounds on the
network activations computed using Bunel et al. (2020). For this
specification, we consider a decomposition of the network such that the
intermediate inner problems can be solved in closed form with linear
multipliers (See Appendix B.6). We leverage Proposition 2 to solve the inner
problem corresponding to the specification. For this specification, with
linear multipliers our approach takes on average under 6 minutes to converge
(per problem instance). See Appendix D for details.

##### Results:

For verification-agnostic networks trained using ACET (Hein et al., 2019) (in
contrast to Bitterwolf et al. (2020a), where they are trained to be
verifiable), we are able to provide stronger guarantees that are closer to the
adversarial bound across settings (See Table 3).

###  6.3 Robust OOD Detection in Stochastic NNs

We now demonstrate results on stochastic neural networks, where the model
parameters are random variables.

##### Verification Task:

We focus on two special cases of practical interest. First, we consider
Bayesian neural networks (BNNs) whose posterior distribution is a truncated
Gaussian distribution. In practice, we re-use an MLP with two hidden layers of
128 units from (Wicker et al., 2020) and truncate their Gaussian posterior
distribution to three standard deviations. As noted earlier, this truncation
can be avoided if guarantees that hold in probability only are acceptable. We
refer to this model as Gaussian-MLP. Second, we consider MC-Dropout networks
with the LeNet architecture (as in Gal & Ghahramani (2016)). We refer to this
model as Dropout-CNN.

##### Method:

In this set of experiments, we use linear functional multipliers. We use the
previously described LP-BP bound-propagation for inner problems where there is
no bound on the parameters. For inner problems that do depend on parameter
bounds, LP-BP is no longer applicable (because of the bilinearity) and we use
instead the bilinear IBP method described in Wicker et al. (2020). We will
refer to this “mixed” bound-propagation method as BP. The corresponding bounds
are obtained with probability 1, and if we were to relax these guarantees to
only hold up to some probability lower than 1, we note that the method of
(Wicker et al., 2020) would offer tighter bounds. In addition, we note that
the “LBP” bound-propagation method from (Wicker et al., 2020) could also be
applied in this experiment, and it would likely benefit both the baseline and
our method. The specification is the same as in the previous subsection: the
models are originally trained on MNIST, and we wish to verify robust OOD
detection on EMNIST. The intermediate inner problems are solved in a manner
similar to the previous subsection. For the final problem corresponding to the
specification, we use projected gradient ascent during training during
optimization and use the exact solver from Appendix B.2.2 to solve the outer
layer problem during evaluation. This experiment is run over 500 samples from
EMNIST.

##### Results:

We note that the running time of both bound-propagation and the functional
Lagrangian method are similar: bound-propagation takes 3s to 30s according to
the setting, while the functional Lagrangian method runs in about 30s in all
settings. The results are presented in Table 4: the functional Lagrangian
framework predominantly outperforms the bound-propagation based baseline.

Model | Perturbation | GAUC (Lower Bound) | AAUC
---|---|---|---
Radius | BP | Ours | Upper Bound
Gaussian-MLP | 0.01 | 0.59 | 0.68 | 0.87
0.03 | 0.46 | 0.54 | 0.89
0.05 | 0.28 | 0.37 | 0.89
Dropout-CNN | 0.01 | 0.44 | 0.28 | 0.72
0.03 | 0.016 | 0.12 | 0.58
0.05 | 0.0 | 0.022 | 0.43

Table 4: Guaranteed Area Under the Curve (GAUC) and Adversarial Area Under the
Curve (AAUC) values for robust OOD detection with stochastic networks. The
linear functional multipliers yield higher GAUC than BP in all settings but
one, where $\epsilon$ has a very low value.

###  6.4 Distributionally Robust OOD Detection

##### Verification Task:

We bound the largest softmax probability across all labels for OOD inputs over
noisy perturbations of the input where the noise is drawn from a family of
distributions with only two constraints for each $p\in\mathcal{P}_{noise}$:
$\omega\in[-\epsilon,\epsilon]\text{ with prob.
}1.,\displaystyle\operatorname*{\mathbb{E}}_{\omega\sim
p}\left[{\exp\left({\omega
t}\right)}\right]\leq\exp\left({\frac{\sigma^{2}t^{2}}{2}}\right)$, for given
constraints $\epsilon,\sigma>0$. The first constraint corresponds to a
restriction on the support of the noise distribution and the second constraint
corresponds to requiring that the noise distribution is sub-Gaussian with
parameter $\sigma$. We note that no prior verification methods, to the best of
our knowledge, address this setting. Thus, as a baseline, we use verification
methods that only use the support of the distribution and perform verification
with respect to worst-case noise realizations within these bounds.

##### Method:

We use networks trained in the same manner as section 6.2. We use the our
framework with multipliers $\lambda_{k}\left({x}\right)=\theta_{k}^{T}x$ for
$k>1$ and linear-exponential multipliers for the input layer:
$\lambda_{1}\left({x}\right)=\theta_{1}^{T}x+\exp\left({\gamma_{1}^{T}x+\kappa_{1}}\right).$
As baselines, we consider a deterministic functional Lagrangian setting with
linear multipliers (as in section 6.2) that only uses the bounds on the
support of the noise distribution $(-\epsilon,\epsilon)$ and an LP-BP baseline
that uses the same information. The inner maximization of $g_{k}$ for $K-1\geq
k\geq 2$ can be done in closed form and we use the approach described in
propositions 1 to solve $\max g_{0},\max g_{1}$ and 2 to solve $\max g_{K}$.

##### Results:

The AUC curves are plotted in Figure 1. The stochastic formulation with
linear-exponential multipliers achieves higher guaranteed AUC, showing the
value of accounting for the noise distribution, instead of relying only on
bounds on the input noise. Computational costs here are similar to those from
Section 6.2.

Figure 1: AUC of distributionally robust OOD Detection on MNIST vs EMNIST
($\epsilon=0.04,\sigma=0.1$). Linear refers to the deterministic formulation
with linear multipliers, LP-BP to bound propagation using linear programming
and stochastic refers to the stochastic formulation using the linear
exponential multipliers.

##  7 Conclusions

We have developed a general framework for verifying probabilistic
specifications. We improve upon prior work in the context of probabilistic
specifications and demonstrate that the framework offers a trade-off between
tightness of verification and computational complexity. We believe that
specific instantiations of the framework can be explored in greater depth,
developing scalable optimization techniques to solve (8) better for various
networks and specifications.

## References

  * Andriushchenko et al. (2020) Andriushchenko, M., Croce, F., Flammarion, N., and Hein, M.  Square attack: a query-efficient black-box adversarial attack via random search.  _European Conference on Computer Vision_ , 2020.
  * Bertsekas (2015) Bertsekas, D. P.  _Convex optimization algorithms_.  Athena Scientific Belmont, 2015.
  * Bitterwolf et al. (2020a) Bitterwolf, J., Meinke, A., and Hein, M.  Certifiably adversarially robust detection of out-of-distribution data.  _Advances in Neural Information Processing Systems_ , 33, 2020a.
  * Bitterwolf et al. (2020b) Bitterwolf, J., Meinke, A., and Hein, M.  Provable worst case guarantees for the detection of out-of-distribution data.  _arXiv preprint arXiv:2007.08473_ , 2020b.
  * Bunel et al. (2020) Bunel, R., Hinder, O., Bhojanapalli, S., et al.  An efficient nonconvex reformulation of stagewise convex optimization problems.  _arXiv preprint arXiv:2010.14322_ , 2020.
  * Cheng et al. (2017) Cheng, C.-H., Nührenberg, G., and Ruess, H.  Maximum resilience of artificial neural networks.  _International Symposium on Automated Technology for Verification and Analysis_ , 2017.
  * Cohen et al. (2017) Cohen, G., Afshar, S., Tapson, J., and van Schaik, A.  Emnist: an extension of mnist to handwritten letters, 2017.
  * Dathathri et al. (2020) Dathathri, S., Dvijotham, K., Kurakin, A., Raghunathan, A., Uesato, J., Bunel, R., Shankar, S., Steinhardt, J., Goodfellow, I., Liang, P., et al.  Enabling certification of verification-agnostic networks via memory-efficient semidefinite programming.  _arXiv preprint arXiv:2010.11645_ , 2020.
  * Dvijotham et al. (2018) Dvijotham, K., Stanforth, R., Gowal, S., Mann, T. A., and Kohli, P.  A dual approach to scalable verification of deep networks.  In _UAI_ , volume 1, pp. 2, 2018.
  * Ehlers (2017) Ehlers, R.  Formal verification of piece-wise linear feed-forward neural networks.  In _International Symposium on Automated Technology for Verification and Analysis_ , pp. 269–286. Springer, 2017.
  * Fazlyab et al. (2019) Fazlyab, M., Morari, M., and Pappas, G. J.  Probabilistic verification and reachability analysis of neural networks via semidefinite programming.  In _2019 IEEE 58th Conference on Decision and Control (CDC)_ , pp. 2726–2731. IEEE, 2019.
  * Feizollahi et al. (2017) Feizollahi, M. J., Ahmed, S., and Sun, A.  Exact augmented lagrangian duality for mixed integer linear programming.  _Mathematical Programming_ , 161(1-2):365–387, 2017.
  * Gal & Ghahramani (2016) Gal, Y. and Ghahramani, Z.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In _international conference on machine learning_ , pp. 1050–1059. PMLR, 2016.
  * Gamrath et al. (2020a) Gamrath, G., Anderson, D., Bestuzheva, K., Chen, W.-K., Eifler, L., Gasse, M., Gemander, P., Gleixner, A., Gottwald, L., Halbig, K., Hendel, G., Hojny, C., Koch, T., Le Bodic, P., Maher, S. J., Matter, F., Miltenberger, M., Mühmer, E., Müller, B., Pfetsch, M. E., Schlösser, F., Serrano, F., Shinano, Y., Tawfik, C., Vigerske, S., Wegscheider, F., Weninger, D., and Witzig, J.  The SCIP Optimization Suite 7.0.  Technical report, Optimization Online, March 2020a.  URL http://www.optimization-online.org/DB_HTML/2020/03/7705.html.
  * Gamrath et al. (2020b) Gamrath, G., Anderson, D., Bestuzheva, K., Chen, W.-K., Eifler, L., Gasse, M., Gemander, P., Gleixner, A., Gottwald, L., Halbig, K., Hendel, G., Hojny, C., Koch, T., Le Bodic, P., Maher, S. J., Matter, F., Miltenberger, M., Mühmer, E., Müller, B., Pfetsch, M. E., Schlösser, F., Serrano, F., Shinano, Y., Tawfik, C., Vigerske, S., Wegscheider, F., Weninger, D., and Witzig, J.  The SCIP Optimization Suite 7.0.  ZIB-Report 20-10, Zuse Institute Berlin, March 2020b.  URL http://nbn-resolving.de/urn:nbn:de:0297-zib-78023.
  * Gowal et al. (2019) Gowal, S., Dvijotham, K. D., Stanforth, R., Bunel, R., Qin, C., Uesato, J., Arandjelovic, R., Mann, T., and Kohli, P.  Scalable verified training for provably robust image classification.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pp. 4842–4851, 2019.
  * Grone et al. (1984) Grone, R., Johnson, C. R., Sá, E. M., and Wolkowicz, H.  Positive definite completions of partial hermitian matrices.  _Linear algebra and its applications_ , 58:109–124, 1984\.
  * Hein et al. (2019) Hein, M., Andriushchenko, M., and Bitterwolf, J.  Why relu networks yield high-confidence predictions far away from the training data and how to mitigate the problem.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 41–50, 2019.
  * Katz et al. (2017) Katz, G., Barrett, C., Dill, D. L., Julian, K., and Kochenderfer, M. J.  Reluplex: An efficient smt solver for verifying deep neural networks.  In _International Conference on Computer Aided Verification_ , pp. 97–117. Springer, 2017.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _International Conference on Learning Representations_ , 2015.
  * Krizhevsky (2009) Krizhevsky, A.  Learning multiple layers of features from tiny images.  Technical report, University of Toronto, 2009.
  * LeCun & Cortes (2010) LeCun, Y. and Cortes, C.  MNIST handwritten digit database.  _NIST_ , 2010.  URL http://yann.lecun.com/exdb/mnist/.
  * Liu et al. (2019) Liu, C., Arnon, T., Lazarus, C., Barrett, C., and Kochenderfer, M. J.  Algorithms for verifying deep neural networks.  _arXiv preprint arXiv:1903.06758_ , 2019.
  * Madry et al. (2017) Madry, A., Makelov, A., Schmidt, L., Tsipras, D., and Vladu, A.  Towards deep learning models resistant to adversarial attacks.  _arXiv preprint arXiv:1706.06083_ , 2017.
  * Nedich & Ozdaglar (2008) Nedich, A. and Ozdaglar, A.  A geometric framework for nonconvex optimization duality using augmented lagrangian functions.  _Journal of global optimization_ , 40(4):545–573, 2008.
  * Netzer et al. (2011) Netzer, Y., Wang, T., Coates, A., Bissacco, A., Wu, B., and Ng, A. Y.  Reading digits in natural images with unsupervised feature learning.  In _NIPS Workshop on Deep Learning and Unsupervised Feature Learning 2011_ , 2011.  URL http://ufldl.stanford.edu/housenumbers/nips2011_housenumbers.pdf.
  * Qin et al. (2019) Qin, C., Martens, J., Gowal, S., Krishnan, D., Dvijotham, K., Fawzi, A., De, S., Stanforth, R., and Kohli, P.  Adversarial robustness through local linearization, 2019.
  * Raghunathan et al. (2018) Raghunathan, A., Steinhardt, J., and Liang, P. S.  Semidefinite relaxations for certifying robustness to adversarial examples.  In _Advances in Neural Information Processing Systems_ , pp. 10877–10887, 2018.
  * Singh et al. (2018) Singh, G., Gehr, T., Mirman, M., Püschel, M., and Vechev, M. T.  Fast and effective robustness certification.  _NeurIPS_ , 1(4):6, 2018.
  * Tjeng et al. (2017) Tjeng, V., Xiao, K., and Tedrake, R.  Evaluating robustness of neural networks with mixed integer programming.  _arXiv preprint arXiv:1711.07356_ , 2017.
  * Vandenberghe & Andersen (2015) Vandenberghe, L. and Andersen, M. S.  Chordal graphs and semidefinite optimization.  _Foundations and Trends in Optimization_ , 1(4):241–433, 2015.
  * Weng et al. (2019) Weng, L., Chen, P.-Y., Nguyen, L., Squillante, M., Boopathy, A., Oseledets, I., and Daniel, L.  Proven: Verifying robustness of neural networks with a probabilistic approach.  In _International Conference on Machine Learning_ , pp. 6727–6736. PMLR, 2019.
  * Wicker et al. (2020) Wicker, M., Laurenti, L., Patane, A., and Kwiatkowska, M.  Probabilistic safety for bayesian neural networks.  In _Conference on Uncertainty in Artificial Intelligence_ , pp. 1198–1207. PMLR, 2020.
  * Zhang et al. (2018) Zhang, H., Weng, T.-W., Chen, P.-Y., Hsieh, C.-J., and Daniel, L.  Efficient neural network robustness certification with general activation functions.  In _Advances in neural information processing systems_ , pp. 4939–4948, 2018.

## Appendix

##  Appendix A Functional Lagrangian Theorem

### Proof of Theorem 1

###### Proof.

The optimization problem (2) is equivalent to the following optimization
problem:

$\displaystyle\max_{p_{0},p_{1},\ldots,p_{K}}\quad$
$\displaystyle\boldsymbol{\psi}\left({p_{K}}\right)$ (10a) Subject to
$\displaystyle
p_{k+1}\left({y}\right)=\int\pi_{k}\left({y|x}\right)p_{k}\left({x}\right)dx\quad\forall
y\in\mathcal{X}_{k+1}$ $\displaystyle\qquad\qquad\qquad\qquad k=0,\ldots,K-1$
(10b) $\displaystyle p_{0}\in\mathcal{P}_{0}$ (10c)

where $\boldsymbol{\psi}$ is a functional of the output distribution $p_{K}$.
We refer to the space of probability measures on $\mathcal{X}_{k}$ as
$\mathcal{P}_{k}$ (not that for $k=0$, this may not be the whole space of
probability measures but a constrained set of measures depending on the
specification we would like to verify). The dual of this optimization problem
can be written as follows:

$\displaystyle\max_{p_{0}\in\mathcal{P}_{0},p_{1}\in\mathcal{P}_{1},\ldots}$
$\displaystyle\boldsymbol{\psi}\left({p_{K}}\right)-\sum_{k=0}^{K-1}\int\lambda_{k+1}\left({z}\right)p_{k+1}\left({z}\right)dz$
$\displaystyle+\sum_{k=0}^{K-1}\displaystyle\int_{\mathcal{X}_{k+1},\mathcal{X}_{k}}\lambda_{k+1}\left({z}\right)\pi_{k}\left({z|x}\right)p_{k}\left({x}\right)dxdz$

where we assigned Lagrange multipliers $\lambda_{k+1}\left({y}\right)$ for
every $y\in\mathcal{X}_{k+1}$. The above optimization problem can be
decomposed as follows:

$\displaystyle\boldsymbol{\psi}\left({p_{K}}\right)-\int_{\mathcal{X}_{K}}\lambda_{K}\left({x}\right)p_{K}\left({x}\right)dx$
$\displaystyle+\sum_{k=0}^{K-1}\int_{\mathcal{X}_{k}}\left({\int_{\mathcal{X}_{k+1}}\lambda_{k+1}\left({y}\right)\pi_{k}\left({y|x}\right)dy-\lambda_{k}\left({x}\right)}\right)p_{k}\left({x}\right)dx$

This is a separable problem in each $p_{k}$ and since $p_{k}$ is constrained
to be a probability measure, the optimal choice of $p_{k}$ (for
$k=1,\ldots,K-1$) is a $\delta$ measure with probability $1$ assigned to the
the $x\in\mathcal{X}_{k}$ that maximizes

$\int_{\mathcal{X}_{k+1}}\lambda_{k+1}\left({y}\right)\pi_{k}\left({y|x}\right)dy-\lambda_{k}\left({x}\right)$

The optimization over $p_{0}$ can be rewritten as follows:

$\max_{p_{0}\in\mathcal{P}_{0}}\int_{\mathcal{X}_{0}}\left({\int_{\mathcal{X}_{1}}\lambda_{1}\left({y}\right)\pi_{0}\left({y|x}\right)dy}\right)p_{0}\left({x}\right)dx$

The optimization over $p_{K}$ can be rewritten as follows:

$\boldsymbol{\psi}^{\star}\left({\lambda_{K}}\right)=\max_{p_{K}\in\mathcal{P}_{K}}\boldsymbol{\psi}\left({p_{K}}\right)-\int\lambda_{K}\left({x}\right)p_{K}\left({x}\right)dx$

The overall dual problem can be rewritten as

$\displaystyle\boldsymbol{\psi}^{\star}\left({\lambda_{K}}\right)+\sum_{k=1}^{K-1}\max_{x\in\mathcal{X}_{k}}\left({\int_{\mathcal{X}_{k+1}}\lambda_{k+1}\left({y}\right)\pi_{k}\left({y|x}\right)dy-\lambda_{k}\left({x}\right)}\right)$
$\displaystyle\quad+\max_{p_{0}\in\mathcal{P}_{0}}\int_{\mathcal{X}_{0}}\left({\int_{\mathcal{X}_{1}}\lambda_{1}\left({y}\right)\pi_{0}\left({y|x}\right)dy}\right)p_{0}\left({x}\right)dx$

Writing this in terms of expected values, we obtain $g\left({\lambda}\right)$.

Plugging in $\lambda^{\star}$ into $g\left({\lambda}\right)$, all the terms
cancel except the first term which evaluates to

$\displaystyle\max_{p_{0}\in\mathcal{P}_{0}}\operatorname*{\mathbb{E}}_{x\sim
p_{0}}\left[{\operatorname*{\mathbb{E}}_{y\in\pi_{0}\left({x}\right)}\left[{\lambda_{1}^{\star}\left({x}\right)}\right]}\right]$
$\displaystyle=\max_{p_{0}\in\mathcal{P}_{0}}\operatorname*{\mathbb{E}}_{x\sim
p_{0}}\left[{\operatorname*{\mathbb{E}}_{y\in\pi_{0}\left({x}\right)}\left[{\operatorname*{\mathbb{E}}_{z\sim\pi_{1}\left({x}\right)}\left[{\lambda_{2}^{\star}\left({z}\right)}\right]}\right]}\right]$
$\displaystyle\ldots=\max_{p_{0}\in\mathcal{P}_{0}}\boldsymbol{\psi}\left({p_{0}}\right)$

∎

##  Appendix B Inner Maximization Problems

###  B.1 Computation of Expected Values

Since $\pi_{k}\left({x}\right)$ is typically a distribution that one can
sample from easily (as it is required to perform forward inference through the
neural network), estimating this expectation via sampling is a viable option.
However, in order to turn this into verified bounds on the specification, one
needs to appeal to concentration inequalities and the final guarantees would
only be probabilistically valid. We leave this direction for future work.

Instead, we focus on situations where the expectations can be computed in
closed form. In particular, we consider layers of the form

$\pi_{k}\left({x}\right)=ws\left({x}\right)+b,(w,b)\sim p^{w}_{k}$

where $s$ is an element-wise function like ReLU, sigmoid or tanh and $(w,b)$
represents a fully connected or convolutional layer. We consider a general
form of Lagrange multipliers as a sum of quadratic and exponential terms as
follows:

$\lambda\left({x}\right)=q^{T}x+\frac{1}{2}x^{T}Qx+\kappa\exp\left({\gamma^{T}x}\right)$

Let

$\tilde{s}\left({x}\right)=\begin{pmatrix}1\\\
s\left({x}\right)\end{pmatrix},\tilde{Q}=\begin{pmatrix}0&q^{T}\\\
q&Q\end{pmatrix},\tilde{w}=\begin{pmatrix}b&w\end{pmatrix}$

Then,

$\lambda\left({Ws\left({x}\right)+b}\right)=\frac{1}{2}\left({\tilde{s}\left({x}\right)}\right)^{T}\tilde{w}^{T}\tilde{Q}\tilde{w}\tilde{s}\left({x}\right)+\kappa\exp\left({\gamma^{T}\tilde{w}\tilde{s}\left({x}\right)}\right)$

Taking expected values with respect to $\tilde{W}\sim p^{w}_{k}$, we obtain

$\displaystyle\operatorname*{\mathbb{E}}\left[{\lambda\left({ws\left({x}\right)+b}\right)}\right]=\frac{1}{2}\left({\tilde{s}\left({x}\right)}\right)^{T}\operatorname*{\mathbb{E}}_{\tilde{s}\left({x}\right)}\left[{\tilde{w}^{T}\tilde{Q}\tilde{w}}\right]+\kappa\prod_{i,j}\operatorname*{\mathbb{E}}\left[{\exp\left({\gamma_{i}\tilde{w}_{ij}\tilde{s}\left({x_{j}}\right)}\right)}\right]$
(11)

where we have assumed that each element of $\tilde{W}_{ij}$ is independently
distributed. First term in (11) evaluates to

$\frac{1}{2}\text{Trace}\left({\text{Cov}\left({\tilde{w}\tilde{s}\left({x}\right)}\right)\tilde{Q}}\right)+\frac{1}{2}\left({\operatorname*{\mathbb{E}}\left[{\tilde{w}\tilde{s}\left({x}\right)}\right]}\right)^{\top}\tilde{Q}\left({\operatorname*{\mathbb{E}}\left[{\tilde{w}\tilde{s}\left({x}\right)}\right]}\right)$

and the second evaluates to

$\kappa\prod_{i,j}\text{mgf}_{ij}\left({\gamma_{i}\tilde{s}_{j}\left({x}\right)}\right)$

where $\text{Cov}\left({X}\right)$ refers to the covariance matrix of the
random vector and $\text{mgf}_{ij}$ refers to the moment generating function
of the random variable $\tilde{w}_{ij}$:

$\text{mgf}_{ij}\left({\theta}\right)=\operatorname*{\mathbb{E}}\left[{\exp\left({\tilde{w}_{ij}\theta}\right)}\right]$

The details of this computation for various distributions on $\tilde{w}$
(Gaussian posterior, MC-dropout) are worked out below:
_Diagonal Gaussian posterior:_ If we have a BNN with a Gaussian posterior,
$\tilde{w}\sim\mathcal{N}\left({\mu,\mathrm{diag}\left({\sigma^{2}}\right)}\right)$,
where $\mu,\sigma\in\mathbb{R}^{mn}$, let
$\text{Mat}\left({\mu}\right)\in\mathbb{R}^{m\times n}$ denote a reshaped
version of $\mu$. Then, we have

$\displaystyle\operatorname*{\mathbb{E}}\left[{\lambda\left({ws\left({x}\right)+b}\right)}\right]=$
$\displaystyle\frac{1}{2}\text{Trace}\left({\mathrm{diag}\left({\text{Mat}\left({\sigma^{2}}\right)\tilde{s}\left({x}\right)}\right)\tilde{Q}}\right)+\frac{1}{2}{\left({\text{Mat}\left({\mu}\right)\tilde{s}\left({x}\right)}\right)}^{T}\tilde{Q}\left({\text{Mat}\left({\mu}\right)\tilde{s}\left({x}\right)}\right)$
$\displaystyle\quad+\kappa\prod_{i,j}\exp\left({\text{Mat}\left({\mu}\right)_{ij}\tilde{s}\left({x_{j}}\right)\gamma_{i}+\frac{1}{2}\text{Mat}\left({\sigma^{2}}\right)_{ij}\left({\tilde{s}\left({x_{j}}\right)\gamma_{i}}\right)^{2}}\right)$

_MC dropout:_ If we have a BNN with a Gaussian posterior,
$\tilde{w}=\mu\odot\text{Bernoulli}\left({p_{\text{dropout}}}\right)$, where
$\mu\in\mathbb{R}^{mn}$ denotes the weight in the absence of dropout and
$p_{\text{dropout}}\in\mathbb{R}^{mn}$ denotes the probability of dropout.
Then, we have

$\displaystyle\operatorname*{\mathbb{E}}\left[{\lambda\left({ws\left({x}\right)+b}\right)}\right]=$
$\displaystyle\frac{1}{2}\text{Trace}\left({\mathrm{diag}\left({\text{Mat}\left({\mu\odot
p_{\text{dropout}}\odot\left({1-p_{\text{dropout}}}\right)}\right)\tilde{s}\left({x}\right)}\right)\tilde{Q}}\right)+\frac{1}{2}{\left({\text{Mat}\left({\mu\odot
p_{\text{dropout}}}\right)\tilde{s}\left({x}\right)}\right)}^{T}\tilde{Q}\left({\text{Mat}\left({\mu\odot
p_{\text{dropout}}}\right)\tilde{s}\left({x}\right)}\right)$
$\displaystyle\quad+\kappa\prod_{i,j}\left({\text{Mat}\left({p_{\text{dropout}}}\right)_{ij}\exp\left({\text{Mat}\left({\mu}\right)_{ij}\tilde{s}\left({x_{j}}\right)\gamma_{i}}\right)+1-\text{Mat}\left({p_{\text{dropout}}}\right)_{ij}}\right)$

###  B.2 Output Layer Optimization

We describe three options here - a convex optimization formulation (Section
B.2.1) that bounds the optimal value for linear multipliers, an exhaustive
optimization formulation for linear multipliers (Section B.2.2) that scales
exponentially with the number of classes but computes the optimal value and a
more general gradient based optimization that can handle both linear and
quadratic multipliers.

####  B.2.1 Convex Optimization Formulation

##### Proof of Proposition 2:

###### Proof.

We remind that we tackle the problem:

$\max\limits_{x_{K}\in\mathcal{X}_{K}}c(x)-\lambda_{K}(x),$

where $c(x)=\mu^{\top}\text{softmax}(x)$. For $i\in\\{1,\ldots,N-1\\}$, let:

$m_{i}=\max\limits_{x_{K}\in\mathcal{X}_{K}}c(x)-\lambda_{K}(x)\text{ such
that }t_{i}\leq c(x)\leq t_{i+1}.$

We have that:

$\displaystyle m_{i}$
$\displaystyle=\max\limits_{x_{K}\in\mathcal{X}_{K}}c(x)-\lambda_{K}(x)\text{
such that: }t_{i}\leq c(x)\leq t_{i+1},$
$\displaystyle=\max\limits_{t,x_{K}\in\mathcal{X}_{K}}t-\lambda_{K}(x)\text{
such that: }t_{i}\leq t\leq t_{i+1},\>t\leq c(x),$
$\displaystyle=\max\limits_{t,x_{K}\in\mathcal{X}_{K}}t-\lambda_{K}(x)\text{
such that: }t_{i}\leq t\leq t_{i+1},\>t\leq\mu^{\top}\text{softmax}(x),$
$\displaystyle\leq
t_{i+1}+\max\limits_{x_{K}\in\mathcal{X}_{K}}-\lambda_{K}(x)\text{ such that:
}t_{i}\leq\mu^{\top}\text{softmax}(x).$

Using the assumption that $\mu$ is a one-hot vector and the slight abuse of
notation that $x-\mu^{\top}x$ subtracts the scalar $\mu^{\top}x$ to all
elements of $x$, we get that:
$\mu^{\top}\text{softmax}(x)=\frac{\exp(\mu^{\top}x)}{1^{\top}\exp(x)}=\frac{1}{1^{\top}\exp(x-\mu^{\top}x)}$.
Thus we can write:

$m_{i}\leq t_{i+1}+\max\limits_{x_{K}\in\mathcal{X}_{K}}-\lambda_{K}(x)\text{
such that: }t_{i}\left({1^{\top}\exp(x-\mu^{\top}x)}\right)\leq 1.$

In other words, we have that $m_{i}\leq M_{i}+t_{i+1}$ for all
$i\in\\{1,\ldots,N-1\\}$. Since we clearly have

$\max\limits_{x_{K}\in\mathcal{X}_{K}}c(x)-\lambda_{K}(x)\leq\max\limits_{i\in\\{1,\ldots,N-1\\}}m_{i},$

this concludes the proof. ∎

####  B.2.2 Exact Optimization Formulation

We now describe an alternative formulation that leads to the exact global
optimum rather than an upper bound on it. The approach scales exponentially in
the dimension of $x_{K}$, which is the same as the number of classes in the
classification problem being solved. Thus, if the number of classes is small
($O(10)$, for example), this approach is practical. We assume throughout this
section that $c$ is a one-hot vector, i.e, $\mu_{i}=1,\mu_{j}=0\quad\forall
j\neq i$, so that the final layer problem is of the form

$\frac{\exp\left({x_{i}}\right)}{\sum_{j}\exp\left({x_{j}}\right)}-\lambda_{K}\left({x}\right)$

We further assume that $\lambda_{K}$ is a linear function, so that the final
optimization problem can be written as

$\max_{l\leq x\leq
u}\frac{\exp\left({x_{i}}\right)}{\sum_{j}\exp\left({x_{j}}\right)}+\sum_{i}\lambda_{i}x_{i}$

We demonstrate how to solve this problem exactly in Algorithm 2. The Algorithm
performs a loop over $3^{n}$ choices where each variable is allowed to be
fixed to its lower bound, upper bound or allowed to be free between in this
range. For each of the $3^{n}$ choices, the first order optima are attained at
points where the gradient of the objective function with respect to the non-
fixed variables is set to $0$. We note that the zero derivative equations can
be solved efficiently in propositions 3 (when $x_{i}$ is not a fixed variable)
and 4 (when $x_{i}$ is a fixed variable). In both cases, there are at most 2
solutions to the $0$ derivative equations and we evaluate whether they satisfy
the bound constraints and if so, choose the point with higher objective value.
Taking the maximum over all $3^{n}$ choices gives the global optimum.

######  Proposition 3.

For any $C\geq 0$, define the function

$f\left({x}\right)=\frac{\exp\left({x_{i}}\right)}{\sum_{j}\exp\left({x_{j}}\right)+C}+\sum_{i}\lambda_{i}x_{i}$

If $\lambda_{i}\not\in[-.25,0]$ or $\lambda_{j}\leq 0$ for some $j\neq i$, $f$
has no stationary points. If $\lambda_{i}\in[-.25,0]$ and $\lambda_{j}\geq 0$
for all $j\neq i$, there are at most two stationary points and they are as
follows: Define

$\displaystyle
p_{i}=\frac{1+\sqrt{1+4\lambda_{i}}}{2},p_{j}=\frac{2\lambda_{j}}{1+\sqrt{1+4\lambda_{i}}}\quad\forall
j\neq i\quad
p_{i}^{\prime}=\frac{1-\sqrt{1+4\lambda_{i}}}{2},p_{j}^{\prime}=\frac{2\lambda_{j}}{1-\sqrt{1+4\lambda_{i}}}\quad\forall
j\neq i$

If $C>0$:
If $\sum_{i}p_{i}<1$,
$x=\log\left({p}\right)+\log\left({\frac{C}{1-\sum_{i}p_{i}}}\right)$ is a
stationary point of $f$. If $\sum_{i}p_{i}^{\prime}<1$,
$x^{\prime}=\log\left({p^{\prime}}\right)+\log\left({\frac{C}{1-\sum_{i}p^{\prime}_{i}}}\right)$
is a stationary point of $f$.
If $C=0$:
$x=\log\left({p}\right)$ is a stationary point if $\sum_{i}p_{i}=1$ and
$x=\log\left({p^{\prime}}\right)$ is a stationary point if
$\sum_{i}p^{\prime}_{i}=1$.

###### Proof.

Let
$\zeta_{k}=\frac{\exp\left({x_{k}}\right)}{\sum_{j}\exp\left({x_{j}}\right)+C}$.
Then,

$\displaystyle\frac{\partial f}{\partial
x_{k}}=\lambda_{k}-\zeta_{k}\zeta_{i}\quad\forall k\neq i$
$\displaystyle\frac{\partial f}{\partial
x_{k}}=\lambda_{i}+\zeta_{i}-\zeta_{i}^{2}$

Thus, at a stationary point, we must have
$\zeta_{i}^{2}-\zeta_{i}+\lambda_{i}=0$ and
$\zeta_{k}=\frac{\lambda_{k}}{\zeta_{i}}\quad\forall k\neq i$. The choices
$\zeta=p,\zeta=p^{\prime}$ constitute the only two solutions to these
equations, which exist when $\lambda_{i}\in[-.25,0]$. To recover $x$ from
$\zeta$, we use the fact that

$x_{i}=\log\left({\zeta_{i}}\right)+t$

where $t$ satisfies $t\left({\sum_{i}\zeta_{i}}\right)=t-C$ so that
$t=\frac{C}{1-\sum_{i}\zeta_{i}}$. The conditions in the theorem on
$p,p^{\prime}$ ensure that $t>0$ so that $x$ is real valued and a valid
stationary point. If $C=0$, $\zeta_{k}$ have to satisfy $\sum_{k}\zeta_{k}=1$
(from the definition of $\zeta_{k}$) and hence the conditions on
$p,p^{\prime}$ summing to $1$ are required in this case. ∎

Algorithm 2 Solving softmax layer problem via exhaustive enumeration

Inputs: $\lambda\in\mathbb{R}^{n},i^{\star}\in\\{1,2,\ldots,n\\}$

$f_{opt}\left({x}\right)\leftarrow\frac{\exp\left({x[i^{\star}]}\right)}{\sum_{j}\exp\left({x[j]}\right)}+\sum_{j}\lambda[j]x[j]$,
$x^{\star}\leftarrow l$, $f^{\star}\leftarrow f_{opt}\left({x^{\star}}\right)$

for $v\in{[\text{Lower},\text{Upper},\text{Interior}]}^{n}$ do

$\text{nonbinding}[i]\leftarrow\left({v[i]==\text{Interior}}\right)$,
$x[i]\leftarrow\begin{cases}l[i]\text{ if }v[i]=\text{Lower}\\\ u[i]\text{ if
}v[i]=\text{Upper}\end{cases}$ for $i=1,\ldots,n$

$C\leftarrow\displaystyle\sum_{i\text{ such that
}\text{nonbinding}[i]==\text{False}}\exp\left({x[i]}\right)$ for
$i=1,\ldots,n$

if $\text{nonbinding}[i^{\star}]$ then

Use proposition 3 to find the stationary point $x^{s}$ of the function
$f\left({x}\right)=\frac{\exp\left({x[i^{\star}]}\right)}{\displaystyle\sum_{j\text{
such that }\text{nonbinding}[j]}\exp\left({x[j]}\right)+C}+\sum_{j\text{ such
that }\text{nonbinding}[j]}\lambda[j]x[j]$ attaining the largest value of
$f\left({x}\right)$ if it exists.

else

$D\leftarrow\exp\left({x[i^{\star}]}\right)$

Use proposition 4 to find the stationary point $x^{s}$ of the function
$f\left({x}\right)=\frac{D}{\displaystyle\sum_{j\text{ such that
}\text{nonbinding}[j]}\exp\left({x[j]}\right)+C}+\sum_{j\text{ such that
}\text{nonbinding}[j]}\lambda[j]x[j]$ attaining the largest value of
$f\left({x}\right)$ if it exists.

end if

if $x^{s}$ exists and $x^{s}[i]\in[l[i],u[i]]\quad\forall
i:\text{nonbinding}[i]$ then

$x[i]\leftarrow x^{s}[i]\quad\forall i:\text{nonbinding}[i]$.

if $f_{opt}\left({x}\right)>f^{\star}$ then

$x^{\star}\leftarrow x$

$f^{\star}\leftarrow f_{opt}\left({x^{\star}}\right)$

end if

end if

end for

Return $x^{\star},f^{\star}$

######  Proposition 4.

For any $C,D>0$, define the function

$f\left({x}\right)=\frac{D}{\sum_{j}\exp\left({x_{j}}\right)+C}+\sum_{i}\lambda_{i}x_{i}$

$f$ has stationary points if and only if

$\lambda>0,\sum_{i}\lambda_{i}\leq\frac{D}{4C}$

If the above conditions are true, $f$ has exactly two stationary points given
by

$\displaystyle
x_{i}=\log\left({\frac{\lambda_{i}}{D}}\right)+2\log\left({t}\right),t=\frac{1\pm\sqrt{1-\frac{4C}{D}\left({\sum_{i}\lambda_{i}}\right)}}{\frac{2}{D}\left({\sum_{i}\lambda_{i}}\right)}.$

###### Proof.

Let
$\zeta_{k}=\frac{\exp\left({x_{k}}\right)}{{\left({\sum_{j}\exp\left({x_{j}}\right)+C}\right)^{2}}}$.
Then,

$\displaystyle\frac{\partial f}{\partial
x_{k}}=\lambda_{k}-D\zeta_{k}\quad\forall k$

Thus, at a stationary point, we must have $\zeta_{k}=\frac{\lambda_{k}}{D}$ so
that we require $\lambda>0$, otherwise no stationary points exist (since
$\zeta>0$ by definition). To recover $x$ from $\zeta$, we use the fact that

$x_{i}=\log\left({\zeta_{i}}\right)+2\log\left({t}\right)$

where $t$ satisfies $t^{2}\left({\sum_{i}\zeta_{i}}\right)=t-C$ so that
$t=\frac{1\pm\sqrt{1-4C\left({\sum_{i}\zeta_{i}}\right)}}{2\left({\sum_{i}\zeta_{i}}\right)}$.
The conditions in the theorem ensure that the solutions for $t$ are real
valued. ∎

####  B.2.3 Optimization with Quadratic Multipliers

When we have diagonal quadratic multipliers as described in Section 5 under
the paragraph “Outer Layer Problem’. Then, the optimization problem can be
rewritten as:

$\displaystyle\max_{l\leq x\leq
u}\frac{\sum_{i}\mu_{i}\exp\left({x_{i}}\right)}{\sum_{i}\exp\left({x_{i}}\right)}-\sum_{i}\alpha_{i}x_{i}-\beta_{i}x_{i}^{2}$
(12)

The following theorem shows how to bound this above:

######  Proposition 5.

Let $[t_{1},t_{2},\ldots,t_{N}]$ be a sequence of real numbers such that
$t_{1}=\min_{i}\mu_{i}$ and $t_{N}=\max_{i}\mu_{i}$. For $i=1,\ldots,N-1$, let
$M_{i}$ be:

$\displaystyle\min_{\theta}\max_{t\in\\{t_{j},t_{j+1}\\}}t+\sum_{i}\max_{l_{i}\leq
z\leq u_{i}}(\mu_{i}-t)\theta\exp\left({z}\right)-\alpha_{j}z-\beta_{j}z^{2}$

Then $\max_{1\leq i\leq N-1}M_{j}$ is an upper bound on the optimal value (9).

###### Proof.

We begin with the formulation in (12). We have that

$t_{1}=\min_{i}\mu_{i}\leq\frac{\sum_{i}\mu_{i}\exp\left({x_{i}}\right)}{\sum_{i}\exp\left({x_{i}}\right)}\leq\max_{i}\mu_{i}=t_{N}$

Thus, the feasible set can be partitioned into $N-1$ parts with

$t_{j}\leq\frac{\sum_{i}\mu_{i}\exp\left({x_{i}}\right)}{\sum_{i}\exp\left({x_{i}}\right)}\leq
t_{j+1}$

where $j=1,\ldots N-1$. We bound the optimization problem over each partition
as follows

$\displaystyle\max_{\begin{subarray}{c}l\leq x\leq u\\\ t_{j}\leq t\leq
t_{j+1}\end{subarray}}$ $\displaystyle
t-\sum_{i}\alpha_{i}x_{i}-\beta_{i}x_{i}^{2}$ Subject to $\displaystyle
t\left({\sum_{i}\exp\left({x_{i}}\right)}\right)=\sum_{i}\mu_{i}\exp\left({x_{i}}\right)$

For any fixed value of $t$, we can dualize the constraint to obtain that the
maximum over $x$ is bounded above by:

$\displaystyle t+\max_{l\leq x\leq
u}\theta\left({\sum_{i}\left({\mu_{i}-t}\right)\exp\left({x_{i}}\right)}\right)-\sum_{i}\alpha_{i}x_{i}-\beta_{i}x_{i}^{2}$
$\displaystyle=t+\sum_{i}\max_{z\in[l_{i},u_{i}]}\left({\mu_{i}-t}\right)\theta\exp\left({z}\right)-\alpha_{i}z-\beta_{i}z^{2}$

for any $\theta\in\mathbb{R}$. The objective is convex in $t$, so that the
maximum over $t$ is obtained at one of the bounds on $t$ \- $t_{j}$ or
$t_{j+1}$. Thus, we obtain an upper bound of

$\max_{t\in\\{t_{j},t_{j+1}\\}}t+\sum_{i}\max_{z\in[l_{i},u_{i}]}\left({\mu_{i}-t}\right)\theta\exp\left({z}\right)-\alpha_{i}z-\beta_{i}z^{2}$

which is valid for any $\theta\in\mathbb{R}$. Minimizing this over $\theta$,
we obtain

$M_{j}=\min_{\theta}\max_{t\in\\{t_{j},t_{j+1}\\}}t+\sum_{i}\max_{z\in[l_{i},u_{i}]}\left({\mu_{i}-t}\right)\theta\exp\left({z}\right)-\alpha_{i}z-\beta_{i}z^{2}$

Taking the maximum over $j=1,\ldots,N-1$, we obtain a bound on (12). ∎

The scalar maximization problems over $z$ can be solved using the following
observation: The second derivative of the objective is
$\theta\left({\mu_{i}-t}\right)\exp\left({z}\right)-2\beta_{i}$ so that the
objective is convex when
$\theta\left({\mu_{i}-t}\right)\exp\left({z}\right)\geq 2\beta_{i}$ and
concave otherwise. Thus, if
$\frac{2\beta_{i}}{\theta\left({\mu_{i}-t}\right)}>0$ (otherwise the objective
is convex or concave over the entire feasible set and can be optimized
efficiently), we can compute
$z^{\ast}=\log\left({\frac{2\beta_{i}}{\theta\left({\mu_{i}-t}\right)}}\right)$
and split the feasible set based on this value. On the convex part, the
maximum is obtained at one of the bounds. For the concave part, the problem
can be solved using binary search efficiently.

####  B.2.4 Final Layer of Robust OOD Detection for Stochastic NNs

We remind that we use linear multipliers, thus the inner problem can be
formulated as:

$\max_{x_{K}\in\mathcal{X}_{K}}\mu^{\top}\texttt{softmax}(x_{K})+\nu^{\top}x_{K},$
(13)

where $\mu$ is a one-hot encoded vector and $\nu$ is a real-valued vector.

##### Projected Gradient Ascent (Training):

  * •

Hyper-parameters: we use the Adam optimizer (Kingma & Ba, 2015), with a
learning-rate of 1.0 and a maximum of 1000 iterations.

  * •

Stopping criterion: when all coordinates have either zero gradient, or are at
a boundary with the gradient pointing outwards of the feasible set.

  * •

In order to help the gradient method find the global maximum, we use a
heuristic for initialization, which consists of using the following two
starting points for the maximization (and then to take the best of the
corresponding two solutions found):

    1. 1.

Ignore affine part ($\nu=0$), which gives a solution in closed form: set
$x_{K}$ at its upper bound at the coordinate where $\mu$ is 1, and at its
lower bound elsewhere.

    2. 2.

Ignore softmax part ($\mu=0$), which also gives a solution in closed form: set
$x_{K}$ at its upper bound at the coordinates where $\nu\geq 0$, and at its
lower bound elsewhere.

##### Convex Upper Bound (Evaluation):

We use Algorithm 2 at evaluation time, which solves the maximization exactly.

###  B.3 Solving the Input Layer Problems with Linear Exponential Multipliers

We recall the setting from Proposition 1. Let
$\lambda_{1}\left({x}\right)=\alpha^{T}x+\exp\left({\gamma^{T}x+\kappa}\right)$,
$\lambda_{2}\left({x}\right)=\beta^{T}x$,
$g_{0}^{\star}=\max_{p_{0}\in\mathcal{P}_{0}}g_{0}\left({p_{0},\lambda_{1}}\right)$,
and
$g_{1}^{\star}=\max_{x\in\mathcal{X}_{1}}g_{0}\left({x,\lambda_{1},\lambda_{2}}\right)$.

######  Proposition 6.

In the setting described above, and with $s$ as the element-wise activation
function:

$\displaystyle
g_{0}^{\star}\leq\alpha^{T}\left({w\mu+b}\right)+\exp\left({\frac{\left\|{w^{T}\gamma}\right\|^{2}\sigma^{2}}{2}+\gamma^{T}b+\kappa}\right),$
$\displaystyle
g_{1}^{\star}\leq\max_{\begin{subarray}{c}x\in\mathcal{X}_{2}\\\
z=s\left({x}\right)\end{subarray}}\beta^{T}\left({w_{2}z+b_{2}}\right)-\alpha^{T}x-\exp\left({\gamma^{T}x+\kappa}\right).$

The maximization in the second equation can be bounded by solving the
following convex optimization problem:

$\displaystyle\min_{\eta\in\mathbb{R}^{n},\zeta\in\mathbb{R}_{+}}\zeta\left({\log\left({\zeta}\right)-1-\kappa}\right)+\mathbf{1}^{T}\max\left({\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({l_{2}}\right),\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({u_{2}}\right)}\right)+\sum_{i}s^{\star}\left({\alpha_{i}+\zeta\gamma_{i},\eta_{i},l_{2i},u_{2i}}\right)$

where

$s^{\star}\left({a,b,c,d}\right)=\max_{z\in[c,d]}-az-bs\left({z}\right)$

###### Proof.

$\displaystyle\max_{\begin{subarray}{c}x\in\mathcal{X}_{2}\\\
z=s\left({x}\right)\end{subarray}}\beta^{T}\left({w_{2}z+b_{2}}\right)-\alpha^{T}x-\exp\left({\alpha^{T}x+\kappa}\right)$
$\displaystyle\leq\min_{\eta}\max_{x\in\mathcal{X}_{2},z\in
s\left({\mathcal{X}_{2}}\right)}\eta^{T}\left({z-s\left({x}\right)}\right)+\beta^{T}\left({w_{2}z+b_{2}}\right)-\alpha^{T}x-\exp\left({\gamma^{T}x+\kappa}\right)$
$\displaystyle\leq\min_{\eta,\zeta}\max_{x\in\mathcal{X}_{2},t}\mathbf{1}^{T}\max\left({\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({l_{2}}\right),\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({u_{2}}\right)}\right)+\beta^{T}b_{2}-\alpha^{T}x-\eta^{T}s\left({x}\right)-\exp\left({t}\right)+\zeta\left({t-\gamma^{T}x-\kappa}\right)$
$\displaystyle\leq\min_{\eta,\zeta}\zeta\left({\log\left({\zeta}\right)-1-\kappa}\right)+\mathbf{1}^{T}\max\left({\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({l_{2}}\right),\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({u_{2}}\right)}\right)+\max_{l_{2}\leq x\leq
u_{2}}-\left({\alpha+\zeta\gamma}\right)^{T}x-\eta^{T}s\left({x}\right)$
$\displaystyle\leq\min_{\eta,\zeta}\zeta\left({\log\left({\zeta}\right)-1-\kappa}\right)+\mathbf{1}^{T}\max\left({\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({l_{2}}\right),\left({\eta+w_{2}^{T}\beta}\right)\odot
s\left({u_{2}}\right)}\right)$ $\displaystyle\quad+\sum_{i}\max_{l_{2i}\leq
x_{i}\leq
u_{2i}}-\left({\alpha_{i}+\zeta\gamma_{i}}\right)x_{i}-\eta_{i}s\left({x_{i}}\right)$

∎

###  B.4 MLP Multipliers

For the MLP multilipers $\lambda:R^{n}\rightarrow R$, we consider multipliers
of the form:

$W_{2}^{\top}\texttt{ReLU}(W_{1}x+b_{1})+b_{2}+W_{\textrm{affine}}^{\top}x+b_{\mathrm{affine}},$
(14)

where $W_{1}\in\mathbb{R}^{m\times
n},W_{2}\in\mathbb{R}^{m},b_{1}\in\mathbb{R}^{m},b_{2}\in\mathbb{R},W_{\mathrm{affine}}\in\mathbb{R}^{n}$.
This lets us control the non-linearity in the multiplier by varying the the
value of $m$ to be $0\textrm{ (affine)},5,20$ (Table 2) and so on.

##### Projected Gradient Ascent (Training):

Since the problem is highly non-concave, we use the following heuristics to
try and find the global maximum:

  * •

Black-box attack (1st phase): we use the Square adversarial attack
(Andriushchenko et al., 2020), with 10000 iterations, 300 random restarts and
learning-rate of 0.01.

  * •

Fine-tuning (2nd phase): We then choose the best attack from the restarts, and
employ projected gradient ascent, with a learning-rate of 0.01 and 500
iterations to fine-tune further.

We also find that using a convex combination between this method and a method
that upper bounds the maximization problem improves results empirically.
Instead of optimizing for $\theta$ as in (8), we use an approximation:

$\min_{\theta}\sum_{k=0}^{K}\alpha_{k}g^{\mathrm{PGA}}_{k}\left({\lambda_{k},\lambda_{k+1}}\right)+g^{\mathrm{LP}}_{k}(1-\alpha_{k})\left({\lambda_{k},\lambda_{k+1}}\right),$

where $\alpha_{k}\in\mathbb{R}_{+}$ (a hyperparameter set to 0.9 in our
experiments), $g^{\mathrm{PGA}}_{k}\left({\lambda_{k},\lambda_{k+1}}\right)$
is an approximate lower bound on the maximization problem
$\underset{x_{k}\in\mathcal{X}_{k}}{\max}g_{k}\left({x_{k},\lambda_{k},\lambda_{k+1}}\right)$
computed using gradient ascent, and
$g^{\mathrm{LP}}\left({\lambda_{k},\lambda_{k+1}}\right)$ is an upper bound on
the maximization problem
$\underset{x_{k}\in\mathcal{X}_{k}}{\max}g_{k}\left({x_{k},\lambda_{k},\lambda_{k+1}}\right)$
based on the LP relaxation (Ehlers, 2017). We find that such regularizing
makes the problem smoother, and helps gradient ascent succeed (consistent with
prior work on adversarial training where such regularization makes gradient
ascent more reliable (Qin et al., 2019)).

##### Mixed Integer Programming (Evaluation):

Using MLP multipliers resuls in inner maximization problems akin to standard
verification of MLPs. Thus at evaluation time, we use the standard MIP
formulation (Cheng et al., 2017) to encode the optimization problem so that we
get a verified upper bound. We find that using the combination of projected
gradient ascent and LP relaxation (as described above) at training time helps
this MIP evaluation run faster. Similar regularization has helped MIP solvers
scale to larger networks in prior work – IBP regularization in Gowal et al.
(2019), and $\ell_{1}$ regularization on the weights in Cheng et al. (2017).

###  B.5 Solving Inner Maximization with Quadratic Constraints

With
$\pi_{k}\left({x}\right)=w\mathrm{\texttt{ReLU}}\left({x}\right)+b,(w,b)\sim
p^{w}_{k}$ and

$\lambda_{k}\left({x}\right)=\frac{1}{2}x^{T}Q_{k}x+q_{k}^{T}x$

We can compute the expected value in (11) to obtain

$\operatorname*{\mathbb{E}}\left[{\tilde{w}^{T}\tilde{Q}_{k+1}\tilde{w}}\right]=\begin{pmatrix}0&m^{T}\\\
m&M\end{pmatrix}$

Then, the inner maximization problem is of the form

$\displaystyle\max_{x\in[l,u]}\frac{1}{2}\left({\mathrm{\texttt{ReLU}}\left({x}\right)}\right)^{T}M{\mathrm{\texttt{ReLU}}\left({x}\right)}+m^{T}\mathrm{\texttt{ReLU}}\left({x}\right)-\frac{1}{2}x^{T}Q_{k}x-q_{k}^{T}x$

for some $m,M$. We can re-write this problem as:

$\displaystyle\max_{\begin{subarray}{c}x\in[l,u]\end{subarray}}$
$\displaystyle\frac{1}{2}y^{T}My+m^{T}y-\frac{1}{2}x^{T}Q_{k}x-q_{k}^{T}x$
s.t. $\displaystyle y\odot x=y\odot y,$ $\displaystyle y\geq 0,$
$\displaystyle y\geq x,$

where we introduce an auxiliary variable $y$ to replace ReLU(x). The
additional constraints introduced on the optimization problem correspond to
the constraint $y=\textrm{{ReLU}}(x)$ (Raghunathan et al., 2018). Following
the approach of (Dathathri et al., 2020), this quadratically constrained
quadratic program (QCQP) can be relaxed into an SDP, and then rewritten as an
eigenvalue optimization problem with dual variables $\gamma$. We then jointly
optimize the dual variables $\gamma$, along with $\theta$, the parameters of
the functional lagrangian multipliers.

For the case when the deterministic neural network is decomposed into layers
of the form
$\pi_{k}\left({x}\right)=\mathrm{\texttt{ReLU}}\left({wx+b}\right)$, we can
follow an approach similar to that outlined above where we turn the problem
into a quadratically constrained quadratic program, and then relax into an SDP
and then proceed as described above.

###  B.6 Solving Inner Maximization with Linear Multipliers

In its general form, the objective function of the inner maximization problem
can be expressed as:

$g_{k}(x_{k},\lambda_{k},\lambda_{k+1})=\operatorname*{\mathbb{E}}_{y\sim\pi_{k}\left({x_{k}}\right)}\left[{\lambda_{k+1}\left({y}\right)}\right]-\lambda_{k}\left({x_{k}}\right)$
(15)

We assume that the layer is in the form $y=Ws\left({x}\right)+b$, where $W$
and $b$ are random variables and $s$ is an element-wise activation function.
Then we can rewrite $g_{k}(x_{k},\lambda_{k},\lambda_{k+1})$ as:

$g_{k}(x_{k},\lambda_{k},\lambda_{k+1})=\operatorname*{\mathbb{E}}_{W,b}\left[{\lambda_{k+1}\left({W\max\\{x_{k},0\\}+b}\right)}\right]-\lambda_{k}\left({x_{k}}\right)$
(16)

We now use the assumption that $\lambda_{k+1}$ is linear:
$\lambda_{k+1}:y\mapsto\theta_{k+1}^{\top}y$. Then the problem can
equivalently be written as:

$\begin{split}g_{k}(x_{k},\lambda_{k},\lambda_{k+1})&=\operatorname*{\mathbb{E}}_{W,b}\left[{\theta_{k+1}^{\top}\left({Ws\left({x_{k}}\right)+b}\right)}\right]-\theta_{k}^{T}x_{k}\\\
&=\left({\operatorname*{\mathbb{E}}\left[{W}\right]^{\top}\theta_{k+1}}\right)^{\top}s\left({x_{k}}\right)+\theta_{k+1}^{\top}\operatorname*{\mathbb{E}}\left[{b}\right]-\theta_{k}^{T}x_{k}\\\
&=\theta_{k+1}^{T}\operatorname*{\mathbb{E}}\left[{b}\right]+\sum_{i}\left({\operatorname*{\mathbb{E}}\left[{W}\right]^{\top}\theta_{k+1}}\right)_{i}s\left({x_{i}}\right)-\left({\theta_{k}}\right)_{i}x_{i}\end{split}$
(17)

Maximizing the RHS subject to $l\leq x\leq u$, we obtain

$\theta_{k+1}^{T}\operatorname*{\mathbb{E}}\left[{b}\right]+\sum_{i}\max_{z\in[l_{i},u_{i}]}\left({\operatorname*{\mathbb{E}}\left[{W}\right]^{\top}\theta_{k+1}}\right)_{i}s\left({z}\right)-\left({\theta_{k}}\right)_{i}z$

where the maximization over $z$ can be solved in closed form for most common
activation functions $s$ as shown in Dvijotham et al. (2018).

So we can simply apply the deterministic algorithm to compute the closed-form
solution of this problem.

##  Appendix C Additional Details on Adversarial Robustness Experiments

###  C.1 Cascading

As described in Section B.4, during training we use projected gradient ascent
to optimize the parameters of the functional lagrangian multipliers, and
during evaluation formulate the inner maximization problems as MIPs and use
the SCIP solver (Gamrath et al. (2020b)) to solve the MIPs. However, these
problems are NP-hard and are computationally expensive as noted in Katz et al.
(2017). The MIPs get progressively harder as we increase the complexity of the
multipliers (from having 0 to 20 $\texttt{ReLU}s$ as each ReLU corresponds to
introducing a binary integer variable). Since, not all MIP problem instances
terminate within the time-out period, we consider a cascaded approach to
reporting measurements in Table 2. Recall that the 0 ReLU units case (with the
network decomposed into layers of the form $y(x)=W\max\\{x,0\\}+b$, except the
first layer) corresponds to the linear relaxation from Ehlers (2017)
(Proposition 7), and we have these bounds available to us at the outset. If
the MIP solver times out, we set the bound for the corresponding multiplier to
be $+\infty$. Let $g_{l}$ be the bound from the MIP solver for a multiplier
with $l$ units, and $\bar{g}_{l}$ be the cascaded bound based on which the
results are reported in Table 2. We calculate the bounds and accuracies as
follows:

$\displaystyle\bar{g}_{0}$ $\displaystyle=\max\\{g_{\mathrm{LP}},{g_{0}}\\}$
$\displaystyle\bar{g}_{5}$ $\displaystyle=\max\\{g_{1},g_{5},\bar{g}_{0}\\}$
(18) $\displaystyle\bar{g}_{20}$ $\displaystyle=\max\\{g_{20},\bar{g}_{5}\\}$
(19)

Here we decompose the problem into layers of the form $y(x)=\max\\{Wx+b,0\\}$
(except the final layer) for all multipliers. We use code from Dathathri et
al. (2020) (which do not consider the constraints from the LP relaxation) to
solve the resulting SDPs for quadratic multipliers. The constraints from the
LP relaxation can also be incorporated into the formulation from Dathathri et
al. (2020), by considering new penalty terms corresponding to the LP
relaxation constraints, and carefully setting the dual variables (we do not do
this). We consider:

$\displaystyle\bar{g}_{\mathrm{quadratic}}=\max\\{g_{\mathrm{quadratic}},{g_{\mathrm{LP}}}\\}.$
(20)

We define
$\bar{g}_{\mathrm{best}}=\max\\{\bar{g}_{\mathrm{quadratic}},\bar{g}_{\mathrm{20}}\\}$.

While improvements are guaranteed because of the cascaded approach, we observe
that the empirical improvements are not merely an artifact of the cascading
described above. In Figure 2, we compare $g_{0}$, $g_{20}$ and
$g_{\mathrm{quadratic}}$ on 50 random problem instances where the SCIP solver
terminates (for all methods that rely on SCIP). In this setting, we see that
the non-cascaded bounds show significant improvements when increasing the
multiplier complexity, and moving from $g_{\mathrm{LP}}$ (equivalent to affine
multipliers) to $g_{20}$ and $g_{\mathrm{quadratics}}$.

###  C.2 Curriculum for MLP Functional Multipliers

For multipliers of the form described in (13), the optimization of the
parameters of the Lagrangrian multiplier is nonconvex when $m>0$. We found
that in practice, the optimization is able to attain better bounds with the
following heuristic: i) we first optimize $W_{\mathrm{affine}}$ and
$b_{\mathrm{affine}}$ with the parameters $W_{1},b_{1},W_{2},b_{2}$ set to
zero – this implies that the multiplier at this point is simply a affine
functional. ii) Subsequently, we freeze $W_{\mathrm{affine}}$ and
$b_{\mathrm{affine}}$ at the optimized values, and optimize
$W_{1},b_{1},W_{2},b_{2}$ further. We find that this leads to better bounds
when compared to simultaneously optimizing for all the parameters in the MLP.

###  C.3 Methodology and hyperpameters

##### Inner Optimization

  * •

For quadratic multipliers, we follow the approach outlined in Appendix B.5 to
solve the inner maximization, where we simultaneously optimize both the
Lagrangian parameters corresponding to the first-order SDP approach outlined
in Dathathri et al. (2020), and the parameters of the quadratic functional
multiplier.

  * •

For MLP functionals, we use the approach outlined in Appendix B.4 during
training with a combination of losses corresponding to projected gradient
ascient and the bound on the inner maximixation problem based on the LP
relaxation (Ehlers, 2017). The network PGD-NN has 3 hidden layers, and this
corresponds to the problem being decomposed into four inner maximization
problems. We find that only the problems corresponding to the intermediate
layers need regularization with the LP bound to make projected gradient ascent
more reliable. For evaluation, we use a MIP solver as described in Appendix
B.4.

##### Outer Optimization

  * •

Quadratic Multipliers: Both the inner dual parameters and functional
lagrangian multiplier parameters from the quadratic functionals are
simultaneously optimized using the rmsprop optimizer, with an initial learning
rate of 0.0003, and divided by 10 every 30000 steps. We run the optimizer for
a total of 150000 steps.

  * •

MLP Multipliers: For learning the parameters corresponding to the affine part
of the MLP multipliers ($W_{\textrm{affine}},b_{\textrm{affine}})$, we use the
rmsprop optimizer with an initial learning rate of 0.001 and divided by 10
every 2000 steps. We run the optimizer for a total of 6000 steps. With these
parameters frozen, we then optimize $W_{1},b_{1},W_{2},b_{2}$ with a similar
procedure.

Figure 2: We compare the certified upper bound from an LP relaxation with
those from our approach using i) an MLP with 20 ReLU units, and ii) quadratic
functions as multipliers, for 50 verification instances where the SCIP solver
terminates for the MLP multipliers. Each verification instance corresponds to
an adversarial objective for a given data point and a given target label. We
see that with expressive multipliers, we are able to mostly improve on the
bounds obtained from the LP relaxation. For a few instances, due to the
difficulty in optimization, we observe bounds with the expressive multipliers
that are slightly worse than the LP bound.

##  Appendix D Additional details on Deterministic OOD experiments

(a) MNIST vs EMNIST (OOD), $\epsilon=.01$

(b) SVHN vs CIFAR(OOD), $\epsilon=0.01$

(c) CIFAR vs SVHN (OOD), $\epsilon=\frac{1}{255}$

Figure 3: Comparison of the AUC curves for the LP-BP baseline (guaranteed AUC)
and our approach (guaranteed AUC) with the upper bounding AUC curve based on
PGD (adversarial AUC) evaluated on networks trained using CEDA (Hein et al.,
2019). The AUC score corresponds to the ability of networks trained on i)
MNIST, ii) SVHN and iii) CIFAR-10 to detect OOD data from i) EMNIST (English
letters), ii) CIFAR-10 and iii) SVHN respectively.

##### Inner Optimization:

All inner problems have a closed-form as shown in section B.6, except for the
last one, which is handled as described in section B.2.2.

##### Outer Optimization:

We use the Adam optimizer, with an initial learning rate of 0.001 for the
network trained on MNIST, and 0.0001 for the networks trained on SVHN and
CIFAR-10. We anneal the learning rate by a factor of 0.1 after 40000 and 70000
steps of optimization, and we run the optimization for a total of 100000
steps.

##### CNN architecture and training

We train networks on MNIST, CIFAR-10 and SVHN using the code from
https://gitlab.com/Bitterwolf/GOOD with the CEDA method, and with the default
hyperparameters. We train a CNN with ReLU activations the following layers:

  * •

Layer 1: Convolutional layer with 16 filters of size $4\times 4$, with padding
dimension of $2$ on both sides, and stride $2$

  * •

Layer 2: Convolutional layer with 32 filters of size $4\times 4$, with padding
dimension of $1$ on both sides, and stride $1$

  * •

Layer 3: Fully connected layer with 100 output neurons

  * •

Logits layer with 10 classes

We use the same architecture for all the datasets.

For the experiments in Section 6.4, we use the network described above trained
on MNIST, and we follow a similar procedure for the outer optimization.

##  Appendix E Additional Details on Stochastic Networks Experiments

##### Inner Optimization:

All inner problems have a closed-form as shown in section B.6, except for the
last one, which is handled as described in section B.2.2.

##### Outer Optimization:

We use the Adam optimizer, with a learning-rate that is initialized at 0.001
and divided by 10 every 250 steps. We run the optimization for a total of 1000
steps.

##### Gaussian-MLP:

We use the ReLU MLP from (Wicker et al., 2020) that consists of 2 hidden
layers of 128 units each. The models are available at
https://github.com/matthewwicker/ProbabilisticSafetyforBNNs.

##### Dropout-CNN:

We use the LeNet5 architecture with dropout applied to the last fully
connected layer with a probability of 0.5. To make the bound-propagation
simpler, we do not use max-pooling layers and instead increase the stride of
convolutions.

Figure 4: Max Softmax Probability per sample on Dropout-CNN ($\epsilon=0.03$),
with the samples ordered by increasing ”linear” upper bound.

Figure 5: AUC of OOD Detection on Gaussian-MLP ($\epsilon=0.03$).

Figure 6: AUC of OOD Detection on Dropout-CNN ($\epsilon=0.03$).

##  Appendix F Relationship to prior work

We establish connections between the functional Lagrangian framework and prior
work on deterministic verification techniques based on Lagrangian relaxations
and SDP relaxations.

###  F.1 Relationship to Dvijotham et al. (2018)

We assume that the network layers are deterministic layers of the form:

$\begin{split}\forall\>k,k\text{ is odd
}\>\pi_{k}\left({x}\right)&=w_{k}x+b_{k},\\\ \forall\>k,k\text{ is even
}\>\pi_{k}\left({x}\right)&=s\left({x}\right),\end{split}$ (21)

where $s$ is an element-wise activation function and that the specification
can be written as:

$\boldsymbol{\psi}\left({x_{K}}\right)=c^{T}x_{K}.$ (22)

######  Proposition 7 (Linear Multipliers).

For a verification problem described by equations (21, 22), the functional
Lagrangian framework with linear functional multipliers
$\lambda_{k}\left({x}\right)=\theta_{k}^{T}x$ is equivalent to the Lagrangian
dual approach from (Dvijotham et al., 2018).

###### Proof.

The final layer problem is

$\max_{x_{K}}c^{T}x_{K}-\theta_{K}^{T}x_{K}=\mathbf{1}^{T}\max\left({\left({c-\theta_{K}}\right)\odot
l_{K},\left({c-\theta_{K}}\right)\odot u_{K}}\right)$

For even layers with $k<K$, the optimization problem is

$\displaystyle\max_{x\in[l_{k},u_{k}]}\theta_{k+1}^{T}\left({w_{k}x+b_{k}}\right)-\theta_{k}^{T}x=\theta_{k+1}^{T}b_{k}+\left({w_{k}^{T}\theta_{k+1}-\theta_{k}}\right)^{T}x$
$\displaystyle=\mathbf{1}^{T}\max\left({\left({w_{k}^{T}\theta_{k+1}-\theta_{k}}\right)\odot
l_{k},\left({w_{k}^{T}\theta_{k+1}-\theta_{k}}\right)\odot
u_{k}}\right)+\theta_{k+1}^{T}b_{k}$

For odd layers with $k<K$, the optimization problem is

$\displaystyle\max_{x\in[l_{k},u_{k}]}\theta_{k+1}^{T}s\left({x}\right)-\theta_{k}^{T}x$
$\displaystyle=\sum_{i}\max_{z\in[l_{ki},u_{ki}]}\left({\theta_{k+1}}\right)_{i}s\left({z}\right)-\left({\theta_{k}}\right)_{i}z$

All these computations precisely match those from Dvijotham et al. (2018),
demonstrating the equivalence. ∎

###  F.2 Relationship to Raghunathan et al. (2018)

We assume that the network layers are deterministic layers of the form:

$\begin{split}\forall\>k\>\pi_{k}\left({x}\right)&=\mathrm{\texttt{ReLU}}\left({w_{k}x+b_{k}}\right)\end{split}$
(23)

where $s$ is an element-wise activation function and that the specification
can be written as:

$\boldsymbol{\psi}\left({x_{K}}\right)=c^{T}x_{K}.$ (24)

######  Proposition 8 (Quadratic Multipliers).

For a verification problem described by equations (23, 24), the optimal value
of the Functional Lagrangian Dual with

$\displaystyle\lambda_{k}\left({x}\right)=q_{k}^{T}x+\frac{1}{2}x^{T}Q_{k}x\quad
k=1,\ldots,K-1$ $\displaystyle\lambda_{K}\left({x}\right)=q_{K}^{T}x$

and when an SDP relaxation is used to upper bound the inner maximization
problems over $g_{k}$, is equal to the dual of the SDP relaxation from
(Raghunathan et al., 2018).

###### Proof.

With quadratic multipliers of the form

$\lambda_{k}\left({x}\right)=q_{k}^{T}x+\frac{1}{2}x^{T}Q_{k}x\quad
k=1,\ldots,K-2$

and

$\lambda_{K}\left({x}\right)=q_{K}^{T}x$

the inner maximization problems for the intermediate layers are of the form

$\max_{\begin{subarray}{c}x\in[l,u]\\\
y=\mathrm{\texttt{ReLU}}\left({wx+b}\right)\end{subarray}}\tilde{q}^{T}y+\frac{1}{2}y^{T}\tilde{Q}y-q^{T}x-\frac{1}{2}x^{T}Qx$

where
$l=l_{k},u=u_{k},\tilde{q}=q_{k+1},\tilde{Q}=Q_{k+1},q=q_{k},Q=Q_{k},w=w_{k},b=b_{k}$.
Let $x\in\mathbb{R}^{n},y\in\mathbb{R}^{m}$ ($n$ dimensional input, $m$
dimensional output of the layer). Further, let
$\tilde{l}=l_{k+1},\tilde{u}=u_{k+1}$.

We can relax the above optimization problem to the following Semidefinite
Program (SDP) (following Raghunathan et al. (2018)):

$\displaystyle\max_{P}$
$\displaystyle\tilde{q}^{T}P[y]+\frac{1}{2}\text{Trace}\left({\tilde{Q}P[yy^{T}]}\right)-q^{T}P[x]-\frac{1}{2}\text{Trace}\left({QP[xx^{T}]}\right)$
(25a) Subject to $\displaystyle
P=\begin{pmatrix}1&\left({P[y]}\right)^{T}&\left({P[x]}\right)^{T}\\\
P[y]&P[yy^{T}]&P[xy^{T}]\\\
P[x]&\left({P[xy^{T}]}\right)^{T}&P[xx^{T}]\end{pmatrix}\in\mathbb{S}^{n+m+1},$
(25b) $\displaystyle P\succeq 0,$ (25c)
$\displaystyle\mathrm{diag}\left({P[xx^{T}]-l\left({P[x]}\right)^{T}-P[x]u^{T}+lu^{T}}\right)\leq
0,$ (25d)
$\displaystyle\mathrm{diag}\left({P[yy^{T}]-\tilde{l}\left({P[y]}\right)^{T}-P[y]\tilde{u}^{T}+\tilde{l}\tilde{u}^{T}}\right)\leq
0,$ (25e) $\displaystyle P[y]\geq 0,P[w]\geq wP[x],$ (25f)
$\displaystyle\mathrm{diag}\left({wP[xy^{T}]}\right)+P[y]\odot
b=\mathrm{diag}\left({P[yy^{T}]}\right).$ (25g)

where the final constraint follows from the observation that:

$y\odot\left({y-wx-b}\right)=0.$

The above optimization problem resembles the formulation of Raghunathan et al.
(2018) except that it only involves two adjacent layers rather than all the
layers at once. Let $\Delta_{k}$ denote the feasible set given the constraints
in the above optimization problem. Then, the formulation of Raghunathan et al.
(2018) can be written as:

$\displaystyle\max$ $\displaystyle\,\,c^{T}y_{K}$ (26a) subject to
$\displaystyle P_{k}\in\Delta_{k}\quad k=0,\ldots,K-1,\,\,l_{K}\leq y_{K}\leq
u_{K},$ (26b) $\displaystyle P_{k+1}[xx^{T}]=P_{k}[yy^{T}]\quad
k=0,\ldots,K-2,$ (26c) $\displaystyle P_{k+1}[x]=P_{k}[y]\quad
k=0,\ldots,K-2,$ (26d) $\displaystyle y_{K}=P_{K-1}[y].$ (26e)

Note that in Raghunathan et al. (2018), a single large $P$ matrix is used
whose block-diagonal sub-blocks are $P_{k}$ and the constraint $P\succeq 0$ is
enforced. Due to the matrix completion theorem for SDPs (Grone et al., 1984;
Vandenberghe & Andersen, 2015), it suffices to ensure postitive
semidefiniteness of the sub-blocks rather than the full $P$ matrix.

Dualizing the last three sets of constraints above with Lagrangian multipliers
$\Theta_{k}\in\mathbb{S}^{n_{k}+n_{k+1}+1},\theta_{k}\in\mathbb{R}^{n_{k}}$
and $\theta_{K}\in\mathbb{R}^{n_{K}}$, we obtain the following optimization
problem:

$\displaystyle\max$
$\displaystyle\,\,c^{T}y_{K}+\theta_{K}^{T}\left({-P_{K-1}[y]+y_{K}}\right)+\sum_{k=0}^{K-2}\text{Trace}\left({\Theta_{k}\left({P_{k+1}[xx^{T}]-P_{k}[yy^{T}]}\right)}\right)+\sum_{k=0}^{K-2}\theta_{k}^{T}\left({P_{k+1}[x]-P_{k}[y]}\right)$
subject to $\displaystyle P_{k}\in\Delta_{k}\quad k=0,\ldots,K-1,$
$\displaystyle l_{K}\leq y_{K}\leq u_{K}.$

The objective decomposes over $P_{k}$ and can be rewritten as:

$\displaystyle\max_{l_{K}\leq y_{K}\leq
u_{K}}\left({c+\theta_{K}}\right)^{T}y_{K}+\sum_{k=0}^{K-1}\max_{P_{k}\in\Delta_{k}}\text{Trace}\left({\Theta_{k-1}P_{k}[xx^{T}]}\right)+\theta_{k-1}^{T}P_{k}[x]-\theta_{k}^{T}P_{k}[y]-\text{Trace}\left({\Theta_{k}P_{k}[yy^{T}]}\right)\quad,$
(27)

with the convention that $\Theta_{K-1}=0,\Theta_{-1}=0,\theta_{-1}=0$. If we
set $Q_{k}=-\Theta_{k-1},q_{k}=-\theta_{k-1}$ for $k=1,\ldots,K$, then the
optimization over $P_{k}$ precisely matches the optimization in (25). Further,
since $\lambda_{K}$ is linear, the final layer optimization simply reduces to:

$\max_{l_{K}\leq x_{K}\leq u_{K}}c^{T}x_{K}-q_{K}^{T}x_{K},$

which matches the first term in (27).

Thus, the functional Lagrangian framework with quadratic multipliers
$\lambda_{k}$ for $k=1,\ldots,K-2$ and a linear multiplier for $\lambda_{K}$
precisely matches the Lagrangian dual of (26) and since (26) is a convex
optimization problem, strong duality guarantees that the optimal values must
coincide.

∎