Maximum Resilience of Artificial Neural Networks
Chih-Hong Cheng, Georg Nührenberg, and Harald Ruess

arXiv:1705.01040v2 [cs.LG] 5 Jul 2017

fortiss - An-Institut Technische Universität München
Guerickestr. 25, 80805 Munich, Germany
{cheng,nuehrenberg,ruess}@fortiss.org
Abstract. The deployment of Artificial Neural Networks (ANNs) in
safety-critical applications poses a number of new verification and certification challenges. In particular, for ANN-enabled self-driving vehicles
it is important to establish properties about the resilience of ANNs to
noisy or even maliciously manipulated sensory input. We are addressing
these challenges by defining resilience properties of ANN-based classifiers
as the maximum amount of input or sensor perturbation which is still
tolerated. This problem of computing maximum perturbation bounds
for ANNs is then reduced to solving mixed integer optimization problems (MIP). A number of MIP encoding heuristics are developed for
drastically reducing MIP-solver runtimes, and using parallelization of
MIP-solvers results in an almost linear speed-up in the number (up to
a certain limit) of computing cores in our experiments. We demonstrate
the effectiveness and scalability of our approach by means of computing maximum resilience bounds for a number of ANN benchmark sets
ranging from typical image recognition scenarios to the autonomous maneuvering of robots.

1

Introduction

The deployment of Artificial Neural Networks (ANNs) in safety-critical applications such as medical image processing or semi-autonomous vehicles poses
a number of new assurance, verification, and certification challenges [2,5]. For
ANN-based end-to-end steering control of self-driving cars, for example, it is
important to know how much noisy or even maliciously manipulated sensory
input is tolerated [14]. Here we are addressing these challenges by establishing
maximum and verified bounds for the resilience of given ANNs on these kinds
of input disturbances.
More precisely, we are defining and computing safe perturbation bounds for
multi-class ANN classifiers. This measure compares the relative ratio-ordering
of multiple, so-called softmax output neurons for capturing scenarios where one
only wants to consider inputs that classify to a certain class with high probability. The problem of finding minimal perturbation bounds is reduced to solving a
corresponding mixed-integer programming (MIP). In particular, the encoding of
some non-linear functions such as ReLU and max-pooling nodes require the introduction of integer variables. These integer constraints are commonly handled

by off-the-shelf MIP-solvers such as CPLEX1 which are based on branch-andbound algorithms. In the MIP reduction, a number of nonlinear expressions are
linearized using a variant of the well-known big-M [10] encoding strategy. We
also define a dataflow analysis [7] for generating relatively small big-M as the
basis for speeding up MIP solving. Other important heuristics in encoding the
MIP problem include the usage of solving several substantially simpler MIP
problems for speeding up the overall generation of satisfying instances by the
solver. Lastly, branch-and-bound is run in parallel on a number of computing
cores.
We demonstrate the effectiveness and scalability of our approach and encoding heuristics by computing maximum perturbation bounds for benchmark sets
such as MNIST [15] and agent games [16]. These cases studies include ANNs for
image recognition and for high-level maneuver decisions for autonomous control
of a robot. Using the heuristic encodings outlined above we experienced a speedup of about two orders of magnitude compared with vanilla MIP encodings.
Moreover, parallelization of branch-and-bound [25] on different computing cores
can yield, up to a certain threshold, linear speed-ups using a high-performance
parallelization framework.
The practical advantages of our approach for validating and qualifying ANNs
for safety-relevant applications are manifold. First, perturbation bounds provide a formal interface between sensor sets and ANNs in that they provide a
maximum tolerable bound on possible sensor errors. These assume-guarantee
interfaces therefore form the basis for decoupling the design of sensor sets from
the design of the classifier itself. Second, our method also computes minimally
perturbed inputs of different classification, which might be included into ANN
training sets for potentially improving classification results. Third, maximum
perturbation bounds are a useful measure of the resilience of an ANN towards
(adversarial) perturbation, and also for objectively comparing different ANNs.
Last, large perturbation bounds are intuitively inversely related with the problem of overfitting, that is poor generalization to new inputs, which is a common
issue with ANNs.
An overview of concrete problems and various approaches to the safety of
machine learning is provided in [2]. We compare our results only with work
that is most closely related to ours. Techniques including the generation of test
cases [17,18,9] or strengthening the resistance of a network with respect to adversarial perturbation [19] are used for validating and improving ANNs. In contrast
to our work, these methods do not actually establish verified properties on the
input-output behavior of ANNs. Formal methods-based approaches for verifying ANNs include abstraction-refinement based approaches [20], bounded model
checking for neural network for control problems [23] and neural network verification using SMT solvers or other specialized solvers [21,13,11]. Instead we rely
on solving MIP problems and parallelization of branch-and-bound algorithms.
In contrast to previous approaches we also go beyond verification and solve optimization problems for ANNs for establishing maximum perturbation bounds.
1

https://www-01.ibm.com/software/commerce/optimization/cplex-optimizer/

input layer

1

output layer

hidden layers

1

(l−1)

x0
Fully
connected

1
(3)

x1
x2

(2)
n1

(3)

n1

x1

ReLU

input-output function

Node structure

Type

=1

(l)

(l)
w
0i

(l)
(l−1) w1i
x1

(l)

(l)
ni

xi

(l)

imi

(l)
w
dl−1 i

(l−1)
xdl−1

(l−1)

(1)

n2

FullyConnected
tan−1
4 inputs
(2x2)
Max-Pooling

xd
Fully connected
ReLU

Fig. 1: An illustration of how a
neural network is defined.

4 inputs
(2x2)

xj1

(l−1)

tan−1

x0

(l)
ni

(l−1)

=

j=0

xi

(l) (l−1)

wji xj

(l)

(l−1)
xi = max(x(l−1)
, xj2 ,
j1
(l−1)

xj3

=1

(l)

(l)
w
0i

(l)
(l−1) w1i
x1

xdl−1

where
Pd(l−1)

(l)

(l−1)

xj2

(l−1)
Max-Pooling xj3
(l−1)
xj4

Fully
connected

(l)

xi = max(0, imi )

(l)

(l)

ni

(l)
w
dl−1 i

(l−1)

, xj4

)

(l)

xi = tan−1 (imi )

xi

where
(l)
imi

=

Pd(l−1)
j=0

(l) (l−1)

wji xj

Fig. 2: Input-output function neurons.

These kinds of problems might also be addressed in SMT-based approaches either by using binary search over SMT or by using SMT solvers that support
optimization such as νZ [4], but it is not clear how well these approaches scale
to complex ANNs. Recent work also targets ReLU [13] or application of a single
image [11,3] (point-wise robustness or computing measures by taking samples).
Our proposed resilience measure for ANNs goes beyond [13,11,3] in that it applies to multi-classification network using the softmax descriptor. Moreover, our
proposed measure is a property of the classification network itself rather than
just a property of a single image (as in [11]) or by only taking samples from the
classifier without guarantee (as in [3]).
The paper is structured as follows. Section 2 reviews the foundations of feedforward ANNs. Section 3 presents an encoding of various neurons in terms of
linear constraints. Section 4 defines our measure for quantifying the resilience
of an ANN, that is, its capability to tolerate random or even adversarial input
perturbations. Section 5 summarizes our MIP encoding heuristics for substantially increasing the performance of the MIP-solver in establishing in minimal
perturbation bounds of ANN. Finally, we present the results of some of our experiments in Section 6, and we describe possible improvements and extensions
in Section 7.

2

Preliminaries

We introduce some basic concepts of feed-forward artificial neural networks
(ANN) [1]. These networks consist of a sequence of layers labeled from l =
0, 1, . . . , L, where 0 is the index of the input layer, L is the output layer, and
all other layers are so-called hidden layers. For the purpose of this paper we
assume that each input is of bounded domain. Superscripts (l) are used to index

layer l-specific variables, but these superscripts may be omitted for input layers.
(l)
Layers l are comprised of nodes ni (so-called neurons), for i = 0, 1, . . . , d(l) ,
(l)
where d is the dimension of the layer l. By convention nodes of index 0 have a
constant output 1; these bias nodes are commonly used for encoding activation
(l−1)
thresholds of neurons. In a feed-forward net, nodes nj
of layer l − 1 are con(l)

(l)

nected with nodes ni in layer l by means of directed edges of weight wji . For
the purpose of this paper we are assuming that all weights in a network have
fixed values, since we do not consider re-learning. Figure 1 illustrates a small
feed-forward network structure with four layers, where each layer comes with a
different type of node functions, which are also main ingredients of convolutional
neural networks. These node functions are specified in Figure 2. The first hidden
(1)
layer of the network in Figure 1 is a fully-connected ReLU layer . Node n2 , for
example, computes the weighted linear sum of all inputs from the previous layer
(1)
as im2 , and outputs the maximum of 0 and this weighted sum. The second
hidden layer is using max-pooling for down-sampling an input representation by
(2)
reducing its dimensionality; node n1 , for example, just outputs the maximum
(3)
of its inputs. Node n1 in the output layer applies the sigmoid-shaped tan−1
on the weighted linear input sum.
fully-connected
Given an input to the network these node
softmax
functions are applied successively from layer 0
(L−1)
x
to L − 1 for computing the corresponding net(L)
(L)
(L−1) 1
x1
n1
n1
work output at layer L. For l = 1 to L we
(l)
(l)
use xi to denote the output value of node ni
(L−1)
(L)
(l)
(l)
(L)
(L−1) x2
n2
x2
and xi (a1 , . . . , ad ) denotes the output value xi
n2
for the input a1 , . . . , ad , sometimes abbreviated
(l)
(L)
by xi (a).
(L)
(L−1)
x3
n3
n3
(L−1)
For the purpose of multi-class classification,
x3
outputs in layer L are often transformed into a
probability distribution by means of the softmax Fig. 3: Topological structure
function
for an output layer with 3
(L−1)
e xi
neurons using softmax.
.
(L−1)
P
xj
e
j=1,...,dL
(L)

In this way, the output xi is interpreted as the probability of the input to
(L−1)
(L−1)
(L−1)
= 2, x3
= 3 of the nodes
be in class i. For the inputs x1
= −1, x2
in Figure 3, for example, the corresponding outputs (0.0132, 0.2654, 0.7214) for
(L)
(L)
(L)
(x1 , x2 , x3 ) sum up to 1.

3

Arithmetic Encoding of Artificial Neural Networks

In a first step, we are encoding the behavior of ANNs in terms of linear arithmetic
constraints. In addition to [13] we are also considering tan−1 , max-pooling and
softmax nodes as commonly found in many ANNs in practice. These encodings
are based on the input-output behavior of every node in the network, and the

main challenge is to handle the non-linearities, which are arising from non-linear
activation functions (e.g., ReLU and tan−1 ), max-pooling and softmax nodes.
Constraints for ReLU and tan−1 nodes as defined
in Figure 2 are separated into, first, an equality
(l)
constraint (1) for the intermediate value imi and,
second, several linear constraints for encoding the
non-linear behavior of these nodes.
X
(l)
(l) (l−1)
imi =
wji xj
(1)

[1, 1]
1
[−1, 2]

0.23
[0, 0.33]
−0.2

ReLU

0.1
[−2, −1]

j=0,...,d(l−1)

We now describe the encoding of the non- Fig. 4: Dataflow analysis
(l)
(l)
(l)
linear functions (xi = max(0, imi ) or xi = for bounding computed
(l)
values in a neural network.
tan−1 (imi )).
Encoding ReLU activation function. The non(l)
(l)
linearity in ReLU constraints xi = max(0, imi ) is handled using the well-known
(l)
big-M method [10], which introduces a binary integer variable bi together with
(l)
(l)
(l)
(l)
(l)
a positive constant Mi such that −Mi ≤ imi and xi ≤ Mi for all possible
(l)
(l)
values of imi and xi . A derivation of the following reduction is listed in the
appendix.
(l)

(l)

Proposition 1 xi = max(0, imi ) iff the constraints (2a) to (4b) hold.
(l)

xi ≥ 0
(l)

(2a)
(l)

xi ≥ imi
(l)

(l)
imi
(l)

(l)

(2b)
(l)

imi − bi Mi
+ (1 −

(l)
(l)
bi )Mi

(l)

(l)

≤0

(3a)

≥0

(3b)

(l)

xi ≤ imi + (1 − bi )Mi
(l)
xi

≤

(l)
(l)
bi M i

(4a)
(4b)

The efficiency of a MIP-solver via big-M encoding heavily depends on the size
(l)
of Mi , because MIP-solvers typically relax binary integer variables to realvalued variables, resulting in a weak LP-relaxation for large big-M s. It is there(l)
fore essential to choose relatively small values for Mi . We apply static analysis [7] based on interval arithmetic for propagating the bounded input values
(l)
through the network, as the basis for generating “good” values for Mi .
(l)

(l)

Max-Pooling. The output xi of a max-pooling node is rewritten as xi =
(l−1)
(l−1)
(l−1)
(l−1)
max(im1 , im2 ), where im1 = max(xj1 , xj2 ) and im2 = max(xj3 , xj4 ).
Encoding the max(x1 , x2 ) function into MIP constraints is accomplished by introducing three binary integer variables to encode y = max(x1 , x2 ) using the
big-M method.

Property-directed encoding of softmax. The exponential function in the definition
of softmax, of course, can not be encoded into a linear MIP constraint. However,
using the proposition below, one confirms that if the property to be analyzed
does not consider the concrete value of output values from neurons but only the
ratio ordering, then (1) it suffices to omit the construction of the output layer,
(L)
(L−1)
and (2) one may rewrite the property by replacing each xi by xi
.
Proposition 2 Given a feed-forward ANN with softmax output layer and a constant α > 0, then for all i, j ∈ {1, . . . , d(L) }:
(L)

(L)

(L−1)

xi1 ≥ α xi2 ⇔ xi1

(L−1)

≥ ln(α) + xi2

.

This equivalence is simply derived by using the definition of softmax, multiplying
by the positive denominator, and by applying the logarithm and the resulting
inequality. The derivation is listed in the appendix.
Encoding tan−1 with error bounds. The handling of non-linearity in tan−1 (im)
is based on results in digital signal processing for piece-wise approximating
tan−1 (im) with quadratic constraints and error bounds. In case −1 ≤ im ≤ 1
the quadratic approximation methods (Eq. (7) of [22]) are used, and tan−1 (im)
is approximated by π4 im + 0.273 im(1 − |im|) with a maximum error smaller
than 0.0038. The absolute value |im| in the formula is removed by encoding
case splits between im ≥ 0 and im < 0 using big-M methods. Otherwise, when
considering the case im > 1 or im < −1, the symmetry condition of tan−1 [24]
1
states that (1) if im > 0 then tan−1 (im) + tan−1 ( im
) = π2 , and (2) if im < 0
π
−1
−1 1
then tan (im) + tan ( im ) = − 2 . This implies that we can create a variable
iminv with a constraint that iminv im = 1, i.e., variable iminv is the inverse of
im. By utilizing the fact that −1 ≤ iminv ≤ 1, the value of tan−1 (iminv ) can be
computed by the formula in (i).
Moreover, case splits are encoded using the big-M method as outlined above.
Since quadratic terms are used, our approach for handling tan−1 nodes requires
solving mixed integer quadratic constraint problem (MIQCP) problems.
Using these approximations for tan−1 (imi ), we obtain lower and upper bounds
for the value of the node variable xi , where the interval between lower and upper
bound is determined by the approximation error of tan−1 . Since the approximation error propagates through the network and using lower and upper bounds
instead of an equality constraint relaxes the problem, our method computes approximations for the measure when it is used for ANNs with tan−1 as activation
function.
Pre-processing based on dataflow analysis. We use interval arithmetic to obtain
relatively small values for big-M , in order to avoid a weak LP-relaxation of the
(l)
(l)
(l)
MIP. Interval bounds for the values of xi are denoted by [Lo(xi ), Up(xi )]. We
are assuming that all input values (at layer l = 0) are bounded, and the output
of bias nodes is restricted by the singleton [1, 1] (the value of the bias is given
by the weight of a bias node). Interval bounds for the values of node outputs

(l)

xi are obtained from the interval bounds of connected nodes from the previous
layers by means of interval arithmetic.
P
(l)
(l)
(l) (l−1)
The output xi of ReLU nodes is defined by imi = j=0,...,d(l−1) wji xj
(l)

(l)

and the ReLU function max(0, imi ). Therefore, interval bounds for xi are
(l)
(l)
computed by first considering the interval bounds Lo(imi ) and Up(imi ), which
(l−1)
are determined by weights of the linear sum and the bounds on xj
. The
(l)

(l)

bounds Lo(imi ) and Up(imi ) are obtained from interval arithmetic as follows:


X
(l)
(l)
(l−1)
(l)
(l−1)
Lo(imi ) =
min wij · Lo(xj
), wij · Up(xj
)
j=0,...,d(l−1)

X

(l)

Up(imi ) =

j=0,...,d(l−1)
(l)



(l)
(l−1)
(l)
(l−1)
max wij · Lo(xj
), wij · Up(xj
) .

(l)

(l)

Given Lo(imi ) and Up(imi ) the bounds on xj are derived using the definition
of ReLU, i.e.,
(l)

(l)

(l)

(l)

[Lo(xi ), Up(xi )] = [max(0, Lo(imi )), max(0, Up(imi ))] .
(l)

(l)

Note that if Lo(xi ) ≥ 0 or Up(xi ) ≤ 0 these bounds suffice to determine
which case of the piece-wise linear ReLU function applies. In this way, the con(l)
straints (2)-(4) maybe dropped and the value of xi is directly encoded using
linear constraints, which reduces the number of binary variables.
(l)
In the case of max-pooling nodes, the output xi is simply the maximum
(l−1)
(l−1)
(l−1)
(l−1)
max(xj1 , xj2 , xj3 , xj4 ) of its four inputs. Therefore, the bounds Lox(l)
i
and Upx(l) on the output are given by the maximum of the lower and uppers
i

bounds of the four inputs respectively. Interval bounds of the outputs for tan−1
are obtained using a polynomial approximation for tan−1 (see below). Finally,
the output of softmax nodes is a probability in [0, 1] which might also be further
refined using interval arithmetic. These bounds on softmax nodes, however, are
not used in our encodings, because of the property-driven encoding of softmax
output layers as described previously.

4

Perturbation Bounds

We define concrete measures for quantifying the resilience of multi-classification
neural networks with softmax output neurons. This measure for resilience is defined over all possible inputs of the network. In particular, our developments do
not depend on probability distributions of training and test data as in previous
work [3]. Maximum resilience of these ANNs is obtained by means of solving
corresponding MIP problems (or MIQCPs in the case of tan−1 activation functions).
We illustrate the underlying principles of maximum resilience using examples
from the MNIST database [15] for digit recognition of input images (see Fig. 5).

Image is likely to be
number 0 or 3 but
unlikely to be 5

Image is
classified
to be 5

smallest
perturbation

Inputs that do not
classify to “5”

A
Inputs that classify to “5”

x5( L ) is largest among x ( L)  x ( L) , x ( L)  x ( L )
0
5
8
5
( L)
( L)

B

{x0 ,..., x9 }

Fig. 5: Finding the smallest possible perturbation for a multi-class classifier to loose
confidence.

Fig. 6: Two images A, B that
both classify to number 5.

Input images in MNIST are of dimension 24 × 24 and are represented as a
vector a1 , . . . , a576 . Input layers of ANN-based multi-digit classifiers for MNIST
therefore consist of 576 input neurons, and the output layer is comprised of 10
(L)
(L)
softmax neurons. Let the output x0 , . . . , x9 at the last layer be the computed
probabilities for an input image to be classified to characters ‘0’ to ‘9’.
To formally define a perturbation, we allow each input ai (i = 1, . . . , d) to
have a small disturbance i , so the input after perturbation is (a1 +1 , . . . , ad +d ).
We sometimes use the concise notation of a +  := (a1 + 1 , . . . , ad + d ) for the
perturbed input. The global value of the perturbation is obtained by taking the
sum of the absolute values of each disturbance i , i.e., |1 | + |2 | + . . . + |d |.
Definition 1 (Maximum Perturbation Bound for m-th classifier)
For a given ANN with d(L) neurons in a softmax output layer and given constants
α ≥ 1 and k ∈ {1, . . . , d(L) − 1}, we define the maximum perturbation bound
for the m-th classifier, denoted by Φm ,2 to be the maximum value such that:
(L)

(L)

For all inputs a = (a1 , . . . , ad ) where xm (a) ≥ α · xj (a) on
all other classes j ∈ {1, . . . , d(L) } \ {m}, we have that for all
perturbations  = (1 , 2 , . . . , d ) where |1 |+|2 |+. . .+|d | < Φm ,
there exist at most k − 1 classes j 0 ∈ {0, 1 . . . , d(L) } such that
(L)
(L)
xm (a + ) ≤ xj 0 (a + ).
Intuitively, the bound Φm guarantees that for all inputs that strongly (defined
by α) classify to class m, if the total amount of perturbation is limited to a value
strictly below Φm , then either (1) the perturbed input can still be classified as m,
or (2) the probability of classifying to m is among the k highest probabilities.
Dually, Φm is the smallest value such that there exists an input that originally
classifies to m, for which the computed probability for class m may not be among
the k highest after being perturbed with value greater than or equal to Φm .
2

For clarity, we usually omit the dependency of Φm from α.

Fig. 5 illustrates an example of an MNIST image being perturbed, where the
neural network considers the perturbed image to be ‘0’ or ‘3’ with at least the
probability of being a ‘5’. The “not among the k highest” property is an indicator
that the confidence of classifying to class m has decreased under perturbation,
as the perturbed input can be interpreted as at least k other classes. In our
experiment evaluations below we used the fixed value k = 2.
Constant α ≥ 1 may be interpreted as indicating the level of confidence of
being classified to a class m. When setting α to 1, the analysis takes all inputs for
which the probability of class m is greater than or equal to the probabilities of
the other classes. Since there might exist an image that has the same probability
for all classes, setting α = 1 may result in a maximum perturbation of zero.
Increasing k helps to avoid this effect, because it requires that at most k − 1
other classes have probabilities greater than or equal to the probility of m. By
picking an α > 1 low-confidence inputs are removed and part (II) of Definition 1
forces the perturbation to be greater than zero. E.g., assume if point B in Fig. 6
is classified to ‘5’ with probability 0.35 and to ‘0’ with probability 0.34, then even
by setting α = 1.1, point B will not be considered in the analysis. By setting
α to 25 one already only considers inputs that classifies to m with probability
higher than 0.95.
Provided that Φm can be computed for each class m (as shown below), one
defines a measure for safe perturbation by taking the minimum of all Φm , and
the measure is computed by computing each Φm independently.
Definition 2 (Perturbation Bound for ANN) For an ANN with L layers
and d(L) softmax neurons in the output layer, a given α ≥ 1, k ∈ {1, . . . , d(L) −1},
and Φm the perturbation bound for the m-th classifier of this ANN from Definition 1, the perturbation bound for ANN is defined as Ξ := min(Φ1 , . . . , ΦdL ).
Based on the dual interpretation above of Definition 1 we are now ready to encode
the problem of finding Φm in terms of the following optimization problem, where
a = (a1 , . . . , ad ) and a +  = (a1 + 1 , . . . , ad + d ).
X

minimize

i=1,...,d

subject to

|i |
(L)

_

x(L)
∀i ∈ {1, . . . , dL } \ m
m (a) ≥ αxi (a)
^
(L)
x(L)
m (a + ) ≤ xi (a + )

I ⊆ {1, . . . , dL } \ m
|I| = k

∀i∈I

and subject to constraints (1)-(4) for ANN encoding.
(5)

Proposition 3 For a given α ≥ 1 and k ∈ {1, . . . , d(L) −1}, the optimal value of
the optimization problem (5) as stated above equals Φm . For ANNs using tan−1
problem (5) yields an under-approximation Φ0m ≤ Φm , because the feasible region
is relaxed due to the approximation of tan−1 .
The first set of conjunctive constraints specifies that the input a = (a1 , . . . , ad )
strongly classifies to m (i.e., satisfies condition I in Def. 1), while the second set
of disjunctive constraints specifies that by feeding the image after perturbation,
the neural network outputs that at least k classes in I are more likely (or equally
likely) than class m (i.e., the second condition in Def. 1 is violated). Therefore, for
inputP
a = (a1 , . . . , ad ) and its associated perturbation  = (1 , . . . , P
d ), we have
that i=1,...,d |i | ≥ Φm . By computing
the
minimum
objective
of
i=1,...,d |i |
P
satisfying the constraints we obtain i=1,...,d |i | = Φm .
We now address the following issues in order to transform optimization problem (5) into a MIP: (1) the objective is not linear due to the introduction of the
absolute value function, (2) the non-linearity of softmax due to the function
(L−1)
P
(L−1)
(L)
xi = e xi
/ j=1,...,dL exj
, and (3) the disjunction in the second set of
constraints.
(i) Transforming objectives. Since the objective |1 | + |2 | . . . , |d | in problem (5)
is not linear, we create new variables abs
in optimization problem (6), where
i
i ∈ {1, . . . , d}, such that every abs
is
greater
than i and −i . Whenever the
i
value is minimized, we have that abs
=
|
|.
i
i
(ii) Removing softmax output layer. Optimization problem (5) contains the in(L)
(L)
equality xm (a1 , . . . , ad ) ≥ αxi (a1 , . . . , ad ). It follows from Proposition 2 that
(L−1)
(L−1)
(a1 , . . . , ad ) ≥ ln(α) + xi
(a1 , . . . , ad ) is
replacing this inequality with xm
sufficient, thereby omitting the exponential function.
(iii) Transforming disjunctive constraints. The disjunctive constraint in problem (5) guarantees at least k classifications with probability equal or higher as m.
We rewrite it by introducing aPbinary variable ci for each class i 6= m. Then we
use (1) an integer constraint i=1,...,d,i6=m ci ≥ k to select k classifications and
(2) the big-M method to enforce that if classification i is selected (i.e., ci = 1),
the probability of classifying to i is higher or equal to the probability of classifying to m.
By applying the transformations (i)-(iii) to the optimization problem (5) we
obtain problem (6), which is a MIP, and it follows from Proposition 3 that
maximum perturbations bounds can be obtained by solving the MIP in (6).
Theorem 1 For a given α ≥ 1 and k ∈ {1, . . . , d(L) − 1}, the optimum of the
MIP in (6) equals Φm for ANNs with ReLU nodes and softmax output layer. For
ANNs using tan−1 it yields an under-approximation.

minimize

Φm :=

subject to

P

abs
i∈{1,...,d} i

(L−1)

x(L−1)
(a) ≥ ln(α) + xi
m
X
ci ≥ k

(a)

∀i ∈ {1, . . . , dL } \ m

i∈{1,...,dL }\m
(L−1)

xi

(a + ) ≥ x(L−1)
(a + ) − M (1 − ci )
m

∀i ∈ {1, . . . , dL } \ m

abs
≥ i
i

∀i ∈ {1, . . . , d}

abs
≥ −i
i

∀i ∈ {1, . . . , d}

ci ∈ {0, 1}

∀i ∈ {1, . . . , dL } \ m

and subject to constraints (1)-(4) for ANN encoding.
(6)

5

Heuristic Problem Encodings

We list some simple but essential heuristics for efficiently solving MIP problems
for the verification of ANNs. Notice that these heuristics are not restricted to
computing the resilience of ANNs, and may well be applicable for other verification tasks involving ANNs.
1. Smaller big-M s by looking back at multiple layers. The dataflow analysis
in Section 3 essentially views neurons at the same layer to be independent.
Here we propose a more fine-grained analysis by considering a fixed number of
(l)
predecessor layers at once. Finding the bound for the output of a neuron xi ,
for example, can be understood as solving a substantially smaller MIP problem
by considering neurons from layer l − 1 and l − 2 when considering two preceding
layers. These MIP problems are independent for each node in these layers and
can therefore be solved in parallel. For each node, we first set the upper bound as
a variable to be maximized in the objective, and trigger the MIP-solver to find
such a value. Relations over integer binary variables can be derived by applying
similar techniques. Notice that these MIPs only generate correct lower and upper
bounds if they can be solved to optimality.
2. Branching priorities. This encoding heuristics uses the given structure of
feed-forward ANNs in that binary integer variables originating from lower layers
are prioritized for branching. Intuitively, variables from the first hidden layer
only depend on the input and it influences all other binary integer variables
corresponding to neurons in deeper layers.

3. Constraint generation from samples and solver initialization. For computing Φm on complex systems via MIP, we use the following three-step process.
ini
First, find an input assignment (aini
1 , . . . , ad ) such that the probability of clas(L) ini
(L) ini
ini
sifying to m is α times larger, i.e., xm (a1 , . . . , aini
d ) ≥ αxj (a1 , . . . , ad ) for
(L)
ini
ini
all j = 1, . . . , d , j 6= m. Finding (a1 , . . . , ad ) is equivalent to solving a
substantially simpler MIP problem without introducing variables 1 , . . . , d and
abs
abs
1 , . . . , d . Second, use Eq. (6) to compute the minimum perturbation by conini
sidering the domain to be size 1, i.e., {(aini
1 , . . . , ad )}. As the domain is restricted
ini
ini
to a single input, all variables a1 , . . . , ad in Eq. (6) are replaced by constants
ini
aini
1 , . . . , ad . This also yields substantially simpler MIP problems, and the computed bound is denoted by Φini
m . Third, and finally, initialize the MIP-solver
by using the computed values from steps 1 and 2, such that the search directly starts
solution
with objective
Φini
m . Also, the constraint
P with a feasible
P
P
ini
ini
−Φm ≤ i=1,...,d i ≤ Φm , as i=1,...,d i ≤ i=1,...,d |i | = Φm ≤ Φini
m , can be
further added to restrict the search space.

6

Implementation and Evaluation

We implemented an experimental platform in C++ for verifying and computing perturbation bounds for neural networks, which is based on IBM CPLEX
Optimization Studio 12.7 (academic version) for MIP solving. We used three different benchmark sets as the basis for our evaluations: (1) MNIST3 for number
characterization, (2) agent games4 , and (3) deeptraffic for simulating highway
overtaking scenarios5 . These benchmarks are denoted by IMNIST , IAgent , and
Ideeptraffic respectively, in the following. For each of the benchmarks we created
neural networks with different numbers of hidden layers and numbers of neurons, which are shown in Tables 1 and 2. All the networks were trained using
ConvNetJS [12].
– Agents in agent games have 9 sensors, each pointing into a different direction
and returning the distances to an apple, poison or a wall, which amounts to
the 27 inputs. Neural networks of various size were trained for an agent that
gets rewarded for eating red things (apples) and gets negative reward when
it eats green things (poison).
– deeptraffic is used as a gamified simulation environment for highway traffic.
The controller is trained based on a grid sensor map, and it outputs high-level
driving decisions to be taken such as switch lane, accelerate or decelerate.
– For MNIST digit recognition [15] has 576 input nodes for the pixels of a
gray-scale image, where we trained three networks with different numbers of
neurons in the hidden layers.
3
4
5

http://cs.stanford.edu/people/karpathy/convnetjs/demo/mnist.html
http://cs.stanford.edu/people/karpathy/convnetjs/demo/rldemo.html
http://selfdrivingcars.mit.edu/deeptrafficjs/

ID instance & output m # inputs; # neuδ
status Time(s)
rons in hidden
M = 104
layers
0 IAgent m=0
27; 300
0.025
inf
1.9
0.05
feas 7.2
1 I2x50
576; 100
0.075
inf
477.8
MNIST m=0
2
0.1
inf
t.o.
3 I2x50
576; 100
0.025
inf
516.8
MNIST m=1
4
0.05
feas 0.5
5 I2x50
576; 100
0.025
inf
0.3
MNIST m=3
6
0.05
inf
303.9
7
0.075 feas 0.3
8 I2x50
576; 100
0.025
inf
0.3
MNIST m=8
9
0.05
inf
146.0
10
0.075 feas 1.1
11 I4x50
576; 200
0.025
inf
464.7
MNIST m=0
12
0.05
inf
t.o.
13 I4x50
576; 200
0.025
inf
t.o.
MNIST m=1
14
0.05
feas t.o.
15 I4x50
576; 200
0.025
inf
t.o.
MNIST m=2
16
0.05 unkown t.o.
17 I4x50
576; 200
0.025 feas 2.7
MNIST m=3
18
0.05
feas 12.5
19 I4x50
576; 200
0.025
inf
t.o.
MNIST m=4
20
0.05 unkown t.o.

Time(s)
dataflow
0.1
26.9
186.8
t.o.
763.9
0.3
0.3
405.1
0.4
0.3
193.5
1.2
489.4
t.o.
t.o.
261.4
t.o.
t.o.
2.7
18.8712
t.o.
t.o.

Time(s)
heuristic
1.+2.
n.a.
n.a.
35.1
2015.9
40.5
328.3
18.7
68.9
151.6
16.5
37.2
185.3
38.08
65.5
128.21
3197.6
54.32
t.o.
45.88
115.1
66.43
t.o.

Table 1: Execution time for verifying perturbation problem over a single input
instance. Time out (t.o.) is set to be 1 hour. Agent games turn out to be quite
simple to solve, therefore no heuristics are being applied (n.a.).
In our experimental validation we focus on efficiency gains of our MIP encodings and parallelization for verifying neural networks, and the computation of
perturbation bound by means of the optimization problem stated in Eq. (6).
Evaluation of MIP Encodings. To understand how dataflow analysis and our
heuristic encodings reduce the overall execution time, we have created synthetic
benchmarks where for each example, we only ask for a given input instance
(e.g., an image) that classifies to m, whether the perturbation bound is below δ.
By restricting ourselves to only verify a single input instance and by not minimizing δ, the problem under verification (local robustness related to an input)
is substantially simpler and is similar to those stated in [13,3]. Table 1 gives
a summary over results being evaluated using Google Computing Engine (16
CPU and 60 GB RAM) by only allowing 12 threads to be used. Compared to a
(l)
naı̈ve approach that sets Mi uniformly to a large constant, applying dataflow
analysis can bring benefits for instances that take a longer time to solve. The
first two heuristics we have implemented are useful for solving some very difficult problems. Admittedly, it can also result in longer solutions times for simpler
instances, but as our ultimate goal is for scalability such an issue is currently
minor. More difficult instances (see I4x50
MNIST in Table 1) could only be solved using
heuristic 1. for preprocessing.

2500

solution time in s

2000

1500

1000

500

0

8

16

24
32
40
48
number of threads

56 60 64

Fig. 7: Execution time vs. the number of threads of five test inputs for IMNIST
with  = 0.01.
Net: # input; # neurons in
hidden layers, output m
IRL : 27;15 m := 0

IRL 27;15 m := 3
Ideeptraffic : 30;70 m
Ideeptraffic : 30;70 m
Ideeptraffic : 45;70 m
Ideeptraffic : 60;70 m

:=
:=
:=
:=

0
3
2
2

α
1.1
1.2
1.5
1.7
1.3
4.022
78.0305
13.5258
2.2704

# of parallelization
12
12
12
12
12
360
360
360
360

Φm

Time (s)

0.1537
0.3006
0.7666
1.2730
0.6904
11.3475
69.9109
7.6226
0.8089

0.4
0.3
0.1
0.1
1.5
421.58
86.40
124.46
2246.8

Table 2: Computation time and results for computing the maximum resilience
Φm .
Effects of Parallelization. For IMNIST we further measured the solution time for
local robustness with  = 0.01 for 10 test inputs using 8, 16, 24, 32 and 64 threads
on machines that have at least as many CPUs as we allow CPLEX to have
threads. The results are shown in Figure 7. It is clearly visible that using more
threads can bring a significant speed-up till 32 cores, especially for instances that
cannot be solved fast with few threads. Interestingly, one can also observe that
for this particular problem (200 neurons in hidden layers), increasing the number
of threads from 32 to 64 does not improve performance (many lines just flatten
from 32 cores onwards). However, for some other problems (e.g., 400 neurons in
hidden layers in hidden layers or computing resilience), the parallelization effect
can last longer to some larger number of threads. We suspect that for problems
that have reached a certain level of simplicity, adding additional parallelization
may not further help.
Computing Φm by solving problem (6). Table 2 shows the result of computing
precise Φm . For simpler problems, we can observe from the first 4 rows of Table 2
that the computed Φm increases, when the value of the parameter α increases.

This is a natural consequence - for inputs being classified with higher confidence,
it should allow for more perturbation to bring to ambiguity. Notably, using a
value of α above its maximum makes the problem infeasible, because there does
not exist an input for which the neural network has such high confidence. For
complex problems, by setting α is closer to its maximum (which can be computed
by solving another substantially simpler MIP that maximizes α for all inputs that
classify to class m), one shrinks the complete input space to inputs with high
confidence. Currently, scalability of our approach relies on sometimes setting a
high value of α, as can be observed in the lower part of Table 2.

7

Concluding Remarks

Our definition and computation of maximum perturbation bounds for ANNs using MIP-based optimization is novel. By developing specialized encoding heuristics and using parallelization we demonstrate the scalability and possible applicability of our verification approach for neural networks in real-world applications.
Our verification techniques also allow to formally and quantitatively compare the
resilience of different neural networks. Also, perturbation bounds provide a formal assume-guarantee interface for decoupling the design of sensor sets from the
design of the neural network itself. In our case, the network assumes a maximum
sensor input error for resilience, and the input sensor sets need to be designed to
guarantee the given error bound. These kinds of contract-based interfaces may
form the basis for constructing more modularized safety cases for autonomous
systems.
Nevertheless, we consider the developments in this paper as only a first tiny
step towards realizing the full potential of formal verification techniques for artificial neural networks and their deployment for realizing new safety-critical
functionalities such as self-driving cars. For simplicity we have restricted ourselves to 1-norms for measuring perturbations but other vector norms may, of
course, also be used depending on the specific needs of the application context. Also, the development of specialized MIP solving strategies for verifying
ANNs, which go beyond the encoding heuristics provided in this paper, may
result in considerable efficiency gains. Notice also that the offline verification
approach as presented here is applied a posteriori to fixed and ”fully trained”
networks, whereas real-world networks are usually trained and improved in the
field and during operation. Furthermore, the exact relationship of our perturbation bounds with the common phenomena of over-fitting in a neural network
classifier deserves a closer examination, since perturbation may also be viewed
as generalization from samples. And, finally, investigation of further measures of
the resilience of ANN is needed, as perturbation bounds do not generally cover
the resilience of ANNs to input transformations such as scaling or rotation.

References
1. Y. S. Abu-Mostafa, M. Magdon-Ismail, and H.-T. Lin. Learning from data, volume 4. AMLBook New York, NY, USA:, 2012.

2. D. Amodei, C. Olah, J. Steinhardt, P. Christiano, J. Schulman, and D. Mané.
Concrete problems in ai safety. arXiv preprint arXiv:1606.06565, 2016.
3. O. Bastani, Y. Ioannou, L. Lampropoulos, D. Vytiniotis, A. Nori, and A. Criminisi.
Measuring neural net robustness with constraints. CoRR, abs/1605.07262, 2016.
4. N. Bjørner, A.-D. Phan, and L. Fleckenstein. νZ-An Optimizing SMT Solver. In
TACAS, pages 194–199. Springer, 2015.
5. S. Bhattacharyya, D. Cofer, D. Musliner, J. Mueller, and E. Engstrom. Certification considerations for adaptive systems. In ICUAS, pages 270–279. IEEE, 2015.
6. M. Bojarski, D. D. Testa, D. Dworakowski, B. Firner, B. Flepp, P. Goyal, L. D.
Jackel, M. Monfort, U. Muller, J. Zhang, X. Zhang, J. Zhao, and K. Zieba. End
to end learning for self-driving cars. CoRR, abs/1604.07316, 2016.
7. P. Cousot and R. Cousot. Abstract interpretation: a unified lattice model for static
analysis of programs by construction or approximation of fixpoints. In POPL, pages
238–252. ACM, 1977.
8. G. Dantzig. Linear programming and extensions. Princeton university press, 2016.
9. I. J. Goodfellow, J. Shlens, and C. Szegedy. Explaining and harnessing adversarial
examples. arXiv preprint arXiv:1412.6572, 2014.
10. I. E. Grossmann. Review of nonlinear mixed-integer and disjunctive programming
techniques. Optimization and engineering, 3(3):227–252, 2002.
11. X. Huang, M. Kwiatkowska, S. Wang, and M. Wu. Safety verification of deep
neural networks. CoRR, abs/1610.06940, 2016.
12. A. Karpathy.
Convnetjs: Deep learning in your browser (2014).
URL
http://cs.stanford.edu/people/karpathy/convnetjs, 2014.
13. G. Katz, C. W. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer. Reluplex:
An efficient SMT solver for verifying deep neural networks. CoRR, abs/1702.01135,
2017.
14. A. Kurakin, I. Goodfellow, and S. Bengio. Adversarial examples in the physical
world. arXiv preprint arXiv:1607.02533, 2016.
15. Y. LeCun, C. Cortes, and C. J. Burges. The mnist database of handwritten digits,
1998.
16. V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and
M. Riedmiller. Playing atari with deep reinforcement learning. arXiv preprint
arXiv:1312.5602, 2013.
17. A. Nguyen, J. Yosinski, and J. Clune. Deep neural networks are easily fooled: High
confidence predictions for unrecognizable images. In CPVR, pages 427–436, 2015.
18. N. Papernot, P. McDaniel, I. Goodfellow, S. Jha, Z. B. Celik, and A. Swami. Practical black-box attacks against deep learning systems using adversarial examples.
arXiv preprint arXiv:1602.02697, 2016.
19. N. Papernot, P. McDaniel, X. Wu, S. Jha, and A. Swami. Distillation as a defense
to adversarial perturbations against deep neural networks. In Oakland, pages 582–
597. IEEE, 2016.
20. L. Pulina and A. Tacchella. An abstraction-refinement approach to verification of
artificial neural networks. In CAV, pages 243–257. Springer, 2010.
21. L. Pulina and A. Tacchella. Challenging SMT solvers to verify neural networks.
AI Communications, 25(2):117–135, 2012.
22. S. Rajan, S. Wang, R. Inkol, and A. Joyal. Efficient approximations for the arctangent function. IEEE Signal Processing Magazine, 23(3):108–111, 2006.
23. K. Scheibler, L. Winterer, R. Wimmer and B. Becker Towards Verification of
Artificial Neural Networks. In MBMV, pages 30–40, 2015.

24. A. Ukil, V. H. Shah, and B. Deck. Fast computation of arctangent functions for
embedded applications: A comparative analysis. In ISIE, pages 1206–1211. IEEE,
2011.
25. Y. Xu, T. K. Ralphs, L. Ladányi, and M. J. Saltzman. Computational experience
with a software framework for parallel integer programming. INFORMS Journal
on Computing, 21(3):383–397, 2009.

Appendix
(l)

(l)

Proposition 1 xi = max(0, imi ) iff constraints (2a) to (4b) hold.
First we establish a lemma to assist the proof.
(l)

(l)

Lemma 1 bi = 1 ⇔ imi ≥ 0.
(l)

(l)

Proof. (⇒) Assume bi = 1, then (3a) holds trivially and (3b) implies imi ≥ 0.
(l)
(l)
(⇐) Assume imi ≥ 0, then (3b) holds trivially and (3a) only holds if bi = 1.
Proof. (Prop. 1)
(l)

(l)

First we rewrite the condition xi = max(0, imi ) to allow further processing.
(l)

(l)

xi = max(0, imi )
definition of max

⇐=========⇒
(l)

Replace im

by b

(l)

(l)

= 1 using lemma 1

(l)

i
i
⇐=======
=====
============⇒

(l)

(l)

(l)

(l)

(l)

(l)

(l)

(imi ≥ 0 ⇒ xi = imi ) ∧ (imi < 0 ⇒ xi = 0)

(l)

(l)

(l)

(l)

(bi = 1 ⇒ xi = imi ) ∧ (bi = 0 ⇒ xi = 0)
(l)

(l)

(l)

(⇒) If (bi = 1 ⇒ xi = imi ) ∧ (bi = 0 ⇒ xi = 0) holds, as bi
integer variable, we consider both cases:
(l)

(l)

is a 0 − 1

(l)

(case bi = 1) From the left clause we derive xi = imi . From Lemma 1 we
(l)
(l)
(l)
(l)
(l)
have imi ≥ 0. By injecting bi = 1, xi = imi , and imi ≥ 0 to con(l)
straints (2a) to (4b), all constraints hold due to very large Mi .
(l)
(l)
(case bi = 0) From the right clause we derive xi = 0. From Lemma 1 we have
(l)
(l)
(l)
(l)
imi < 0. By injecting bi = 0, xi = 0, and imi < 0 to constraints (2a)
(l)
to (4b), all constraints hold due to very large Mi .
(⇐) If all constraints in (2a) to (4b) hold, we do case split to consider cases
(l)
(l)
(l)
(l)
(l)
(l)
bi = 0 and bi = 1, and how they make (bi = 1 ⇒ xi = imi ) ∧ (bi = 0 ⇒
(l)
xi = 0) hold.
(l)

(l)

(l)

(case bi = 1) From (1b) and (3a) we know that xi = imi .
(l)
(l)
(case bi = 0) From (1a) and (3b) we know that xi = 0.
(l)

(l)

(l)

(l)

(l)

In both cases, (bi = 1 ⇒ xi = imi ) ∧ (bi = 0 ⇒ xi = 0) holds.

Proposition 2 Given a feed-forward ANN with softmax output layer and a constant α > 0, then for all i, j ∈ {1, . . . , d(L) }:
(L)

(L)

(L−1)

xi1 ≥ α xi2 ⇔ xi1

(L−1)

≥ ln(α) + xi2

.

Proof.
(L)

(L−1)

⇐
⇒
⇐
⇒

P

exi1

j=1,...,dL

(L)

xi1 ≥
(L−1)

exj

(L−1)

xi1

≥
≥

α x i2
(L−1)

αP

exi2

(L−1)

j=1,...,dL

e xj

(L−1)

ln(α) + xi2

