Neural Network Verification

Motivation

For any reasonably complex system, it is necessary to ensure that the system is acting as intended without any unpredictable behaviours that may end up compromising the system and/or the systems surrounding it. There are many safety critical systems whose malfunctioning can cause innumerable casualties and loss of life, thus for these systems we need to reliably prove that they have predictable outputs corresponding to determined inputs.

Formal Verification utilizes a lot of different mathematical tools to demonstrate the “correctness” of some system with respect to some desirable specification or property. The set of tools used for Formal Verification is collectively known as Formal Methods and include such theories like logic, automata, type systems and formal languages. These can be used to check the behaviour of some system against a set of specifications in order to determine correctness of the system. Formal Verification can be imagined as a layer of abstraction on top of systems (both software and hardware) that provides some guarantees as to the functioning of these systems.

Below we will be building up to the application of Formal Methods to the field of Machine Learning in order to Verify Neural Networks.

Note

The following is an effort to condense the Introduction to Neural Network Verification book in order to concisely summarize the material that is relevant to my research.

Introduction

Background

One of the first demonstrations of the possibilities of verification came from Alan Turing’s 1949 paper called ‘Checking a large routine’. In the paper the program for finding the factorial of a number was broken down into a simple sequential set of instructions that could either be either $True$ or $False$. The truth values for each of these assertions were then checked in order to prove that the program is a faithful implementation of the factorial mathematical function. This is known as proving the functional correctness of a program. Although this is the gold standard for demonstrating correctness, it is not possible to do so for Neural Networks as the tasks performed cannot be described mathematically. In a Neural Network we generally test for the following qualities:

• Robustness: This class of properties is violated when small perturbations to the inputs result in changes to the output. It also measures how well a model performs on new and unseen test data.
• Safety: This is a broad class of correctness properties that ensure that the program does not reach a ‘bad’ state. (e.g., a surgical robot should not be operating at 200mph)
• Consistency: This class of properties are violated only when the algorithms are not consistent with the real world. (e.g., tracking an object in free fall should be consistent with the laws of gravitation)

We will mostly be testing for Robustness, as that is the set of properties that are most widely tested for. It can also encapsulate the Safety and Consistency classes (as there is not well-defined boundary between them) if the specifications are designed accordingly.

Note

Alan Turning had also extensively thought about ‘Networks that can Learn’, which can be seen from his 1948 paper Intelligent Machinery where he proposes Boolean NAND networks that can learn over time.

Abstraction of Neural Network

We will not be dealing with the technicalities of the Neural Networks but rather used well-rounded mathematical abstractions that can be rigorously tested over. Since all Neural Networks can be though of as Directed Acyclic Graphs (DAGs), we will be treating them as such. More specifically they will be treated as dataflow graphs of operators over $\mathbb{R}$. The shape of these graphs will determine what specific architecture they belong to. Each node will perform some computation, whose dependencies are the edges. Thus, a neural network will be associated with some function:
$f:{\mathbb{R}}^n\to {\mathbb{R}}^m$

For a Neural Network: $G=(V,E)$, we have:
$V$: finite set of nodes
$E\subseteq V\times V$: set of edges
$V^{in}\subset V$ : input nodes
$V^o\subset V$ : output nodes
$n_v$ : total number of edges whose target is $v$
where, $n=|V^{in}|$ and $m=|V^o|$ for the network.

All neural networks are finite DAGs, there are classes of neural networks called Recurrent Neural Networks that have loops, but the number of iterations will depend on the size of the input. RNNs have self-loops that unroll based on the length of the input, which means that it reliably terminates. So while testing neural networks we will not test for loop termination and the explosion of possible program paths that comes with it.

For our purposes the following conditions must be met in a network:

• all nodes must be reachable from some input node
• every node can reach an output node
• fixed total ordering on edges $E$ and another one on nodes $V$.

Each node $v$ of the neural network is a function of the following form: $f_v:{\mathbb{R}}^{n_v}\to \mathbb{R}$

Where each node takes in some vector, does some computation on it and return a single output number that is passed through an activation function (a non-linear function), for ease of analysis they are treated as belonging to two different nodes. Each of the elements of the vectors are the outputs of previous nodes. These relations can be recursively defined with the base case terminating at the input nodes. Therefore, for every non-input node $v\in V$, we have:

Where each $(v_i,v)$ represents an edge connecting node $v_i$ to $v$. The ordering of these edges and nodes determines the inputs on the basis of which the computations will be done. Each node will have an $out(v)$ function which can be defined as follows: $out(v)=f_v(x_1,...,x_{n_v})$ where $x_i=out(v_i)$ for $i\in \{1,2,...,n_v\}$. This is the recursive definition wherein each input to the node $x_i$ can be defined as: $out(v)=f_v(out(v_1),out(v_2),...,out(v_{n_v}))$
No sequence of operations are defined, only which nodes need what data to perform its computations. A more modified version of these graphs are known as computation graphs. $out(v_i)$ values can be computed in any topological ordering of graph nodes, as it needs to be ensured that all the inputs need to be computed before the target node itself. The total ordering ensures proper reference to the graph nodes and whether computations have been done in their desirable sequence, below the total ordering of nodes helps linearize the nodes into a topologically sorted sequence.

flowchart LR
1(( 1 ))
2(( 2 ))
3(( 3 ))
4(( 4 ))
5(( 5 ))
6(( 6 ))
7(( 7 ))
8(( 8 ))
9(( 9 ))
10(( 10 ))
11(( 11 ))
12(( 12 ))

1 & 2 & 3 & 4 --> 5 & 6 & 7
5 & 6 & 7 --> 8 & 9
8 & 9 --> 10 & 11 & 12

subgraph input
1
2
3
4
end

subgraph hidden1
5
6
7
end

subgraph hidden2
8
9
end

subgraph output
10
11
12
end

All vector computations need to be linear: $f(x)={\sum }_{i=1}^n{c}_i{x}_i+b$ or piece-wise linear:

$$f(x)\{\begin{array}{l}{\sum }_i{c}_i^{[1]}{x}_i+b^{[1]},x\in S_1\\ ⋮\\ {\sum }_i{c}_i^{[m]}{x}_i+b^{[m]},x\in S_m\end{array}$$

where ${\cup }_i{S}_i={\mathbb{R}}^n$ and ${\cap }_i{S}_i=\varnothing$.

Note

Neural Networks are an instance of differential programs.

Defining Specifications

We define a language that specifies some properties about the functioning of a neural net. This will enable us to later on make statements and verify them based on the specifying language.

Specifications are generally of the form:

$$\begin{array}{r}\{precondition\}\\ r\leftarrow f(x)\\ \{postcondition\}\end{array}$$

where both preconditions and post-conditions are statements that specify some property that is adjacent to input and output respectively. Properties dictate the input-output behaviour of the network (and not the internals). Specifications help in quantifying some properties for accurate verification. Each specification can be thought of as being structured in the following way:

$\underset{precondition}{\underbrace{foranyinputsx,y,\dots that\dots }}theneuralnetworkGproduces\underset{pstcondition}{\underbrace{outputthat\dots }}$

Although every possible specification needed for complete verification cannot be made, multiple specifications can be combined together to test for stronger properties. Preconditions are generally predicates or Boolean functions defined over set of variables which act as inputs to the system, and post condition is a Boolean predicate over the variables appearing in precondition ($x_i$) and assigned variables ($r_i$).
For any values of $x_1,\dots ,x_n$ that make the precondition true, let $r_1=f(x_1)$, $r_2=g(x_2)\dots$ , where $f(),g(),\dots$ are the computations on the input, then the post condition must also be True.

$$\begin{array}{r}\{precondition\}\\ r_1\leftarrow f(x_1)\\ r_2\leftarrow g(x_2)\\ ⋮\\ \{postcondition\}\end{array}$$

If the post condition is False, then the correctness property is not true, i.e., the property does not hold.

Consider: c is an actual greyscale image, each element of c is the intensity of a pixel $\in (0,1)$; we can state the following specification about the brightness of c and its corresponding classification:

$$\begin{array}{r}\{|x-c|\le 0.1\}\\ r_1\leftarrow f(x)\\ r_2\leftarrow f(c)\\ \{class(r_1)=class(r_2)\}\end{array}$$

$r_1$ and $r_2$ are vectors whose elements are probabilities for a belonging to a class with labels corresponding to the indexes. $class(r_1)$ and $class(r_2)$ extract the indexes corresponding to the largest elements of the vectors $r_1$ and $r_2$ respectively.

Info

Counterexamples: Valuations of variables in precondition that falsifies the post condition.

flowchart TD

A["Verification"]
B["Constraint Based"]
C["Abstraction Based"]

A --> B & C 

Constraint-Based Verification

Constraint-Based Satisfaction

A correctness property is taken and encoded as a set of constraints, solving which will help us to decide whether the property holds or not.

Let $fv(F)$ be the set of free variables appearing in the formula $F$.

Interpretation: $I$ of $F$ is a map from variables present in $fv(F)$ to either $True$ or $False$. $I(F)$ denotes the formula where the variables have been replaced with their corresponding interpretations.

Example:
Let $F\triangleq (p\wedge q)\neg r$, which means that $F$ is syntactically defined to be equal to the Boolean formula (as opposed to being semantically equivalent).

$$\begin{array}{r}\\ fv(F)=\{p,q,r\}\\ I=\{p\mapsto \text{True, }q\mapsto \text{False, }r\mapsto \text{True}\}\\ I(F)\triangleq (\text{True}\wedge \text{False})\neg \text{True}\end{array}$$

$eval(F)$: denotes the simplest form of $F$ we can get by evaluating repeatedly.

$F$ is satisfiable (SAT) if there exists an $I$ s.t. $eval(I(F))=True$, in which case, $I$ is the model of $F$: $I\models F$.
$I/\models F$ denotes $I$ isn’t a model of $F$.
$I/\models F$ if and only if $I\models \neg F$.
$F$ is unsatisfiable (UNSAT) if $\forall I,eval(I(F))=False$.

Validity: If every interpretation $I$ is a model of $F$, then $F$ is valid.

Boolean satisfiability theories (SAT) can be generalised to more complex theories to include reals, vectors, strings, arrays, etc.; the problem of determining whether statements within these theories are true is known as Satisfiability Modulo Theories or simply SMT. We will be extensively using a first-order logic system called Linear Real Arithmetic (LRA) as it can represent a large class of neural networks and is decidable.
In LRA, each propositional variable is replaced by a linear inequality of the form:

$$\sum _{i=1}^n{c}_i{x}_i+b\le 0\text{ OR }\sum _{i=1}^n{c}_i{x}_i+b<0$$

where $c_i,b\in \mathbb{R}$.

$$(\underset{p}{\underbrace{x+y\le 0}}\wedge \underset{q}{\underbrace{x-2y<10}})\vee \underset{r}{\underbrace{x\ge 100}}$$

An interpretation $I$ of $F$ is an assignment of every free variable to a real number.

Encoding Neural Networks

We need to translate Neural Networks into a formula in LRA such that we can use SMT Solvers (specialised software that tests satisfiability for SMT problems).

$f_v$: Function in node.
$I$: Model.
$F_v$: Encoding for node.
$R_g$: Relational mapping for neural network.

Note

Generally Variables with a subscript $G$ refers to the overall neural network as a graph and those with a subscript v refers to individual nodes or a collection of nodes.

Whole Network can be defined as a binary relation:

$$R_G=\{(a,b)|a\in {\mathbb{R}}^n,b=f_G(a)\}$$

$R_v$, $f_v$ define the same for a single node $v$ in $G$ network.

For a single neuron with corresponding function $f_v:\mathbb{R}\to \mathbb{R}$, both can be defined as:

$$f_v(x)=x+1\text{and}{R}_v=\{(a,a+1)|a\in \mathbb{R}\}$$

Thus the encoding for a node with one input is as follows

$$F_v\triangleq v^o=v^{in,1}+1$$

Models of $F_v$ will be of the form $\{v^{in,1}\mapsto a,v^o\mapsto a+1\}$ and have one-to-one correspondence with the elements of $R_v$.

For two inputs: The formula for the encoding of a node with the function $f(x)=x_1+1.5x_2$ can be represented as

$F_v\triangleq v^o=v^{in,1}+1.5v^{in,2}$

Generalizing encoding for any single node:

Formalizing the operation $f_v$ of some node $v$. We assume that the function $f_v:{\mathbb{R}}^{n_v}\to \mathbb{R}$ is piecewise linear, i.e., of the form

$$f(x)=\{\begin{array}{ll}{\sum }_j{c}_j^1{x}_j+b& \text{ if }{S}_1\\ ⋮& \\ {\sum }_j{c}_j^l{x}_j+b& \text{ if }{S}_l\end{array}$$

where $j$ ranges from $1$ to $n_v$. This generalizes representations to allow any one linear operation corresponding to some predetermined condition $S_i$. Each $S_i$ is defined as a formula in LRA over the elements of $x$. Thus, the encoding for the single node can be written as:

$$F_v\triangleq \underset{i=1}{\overset{l}{\bigwedge }}\left[S_i\Rightarrow \left(v^{\mathrm{o}}=\sum _{j=1}^{n_v}{c}_j^i\cdot v^{\mathrm{in},j}+b^i\right)\right]$$

if statement $S_i$ is $True$, then $v^o$ is equal to the $i^{th}$ equality. The statements are combined together using a conjunction over all possible inputs ($v^{in,j}$) to the node. Each clause joined using the conjunction is in $\text{condition}⟹\text{assignment}$ form providing the functionality of an if statement. So the general way to read this condensed statement is:

$$\text{if }{S}_1\text{ then }\dots \text{ AND if }{S}_2\text{ then }\dots$$

Example (ReLU):

$$relu(x)=\{\begin{array}{lcc}x& \text{if}& x>0\\ o& \text{if}& x\le 0\end{array}$$

The encoding for a simple ReLU function will be:

$$F_v\triangleq (\underset{x>0}{\underbrace{v^{\mathrm{in},1}>0}}\Rightarrow v^0=v^{\mathrm{in},1})\wedge (\underset{x⩽0}{\underbrace{v^{\mathrm{in},1}⩽0}}\Rightarrow v^0=0)$$

In order to ensure that we are making accurate representations of the actual system that we are modelling, we have to ensure that two qualities are always guaranteed for the encodings: Soundness & Completeness.

Soundness - This is to make sure that our model does not miss any behaviour of $f_v$. Let $(\vec{a},b)\in \mathbb{R}$

$$I=\{v^{in,1}\mapsto a_1,\dots ,v^{in,n}\mapsto a_n,v^o\mapsto b\}$$

for any given tuple of $\vec{a}$ and $b$ which is an element of $R_v$, $I$ is a model of $F_v$ : $I\models F_v$. Soundness is the property of only being able to prove “true” things, any analysis about the system that is proven to be true, will be true.

Completeness - Any model of $F_v$ maps to a behaviour of $f_v$, so if for a model of $F_v$:

$$I=\{v^{in,1}\mapsto a_1,\dots ,v^{in,n}\mapsto a_n,v^o\mapsto b\}$$

then $(\vec{a},b)\in \mathbb{R}$. Completeness is the property of being able to prove all true things that can possibly exist in the system. This property can be used to determine counterexamples for any given model.

So far we have seen how to encode a singular node in a larger network. But to encode an entire Neural Net, we need to have the following structure:

flowchart TD
a[Encoding a Neural Network]
b[Encoding a formula for the nodes]
c[Encoding a formula for the edges]
a --> b
a --> c

Encoding the nodes:

for all non-input nodes we have:

$$F_V\triangleq \underset{v\in V/V_{in}}{\bigwedge }{F}_v$$

Since the input nodes do not perform any operations, we are excluding them from the set using $V/V_{in}$.

$$\text{the output of }{v}_1\text{ is }\dots \text{ AND the output of }{v}_2\text{ is }\dots$$

Encoding the edges:

For some node $v\in V/V_{in}$ there exists a total ordering of edges: $(v_1,v),(v_2,v),\dots$ The total ordering tells us which edge feeds into which corresponding input index of the node.

$$F_{o\to v}\triangleq \underset{i=1}{\overset{n}{\bigwedge }}{v}^{in,i}=v^o$$ $$F_E=\underset{v\in V/V_{in}}{\bigwedge }{F}_{o\to v}$$

As per the total ordering of the input edges, $F_{o\to v}$ determines that the output of the $i^{th}$ node from the previous layer will become the $i^{th}$ input of the present node. The $v$‘s present in the below diagram are the sequentially ordered nodes of the previous layer.

More concretely:

$$v_i^{[o]}⟹v^{[in,i]}$$

The entire network can be represented as a conjunction of the edges and nodes:

$$F_G\triangleq F_V\wedge F_E$$

Assuming that we have ordered input nodes in $V_{in}$ Let $(a,b)\in R_G$ and let

$$I=\underset{\text{input}}{\underbrace{\{v_1^o\mapsto a_1,\dots ,v_n^o\mapsto a_n\}}}\cup \underset{\text{output}}{\underbrace{\{v_{n+1}^o\mapsto a_{n+1},\dots ,v_{n+m}^o\mapsto a_{n+m}\}}}$$

Then there exists $I^{\prime }$ such that $I\cup I^{\prime }\models F_G$.

Note

Size of encoding is linear in the size of the Neural Network

Example:

flowchart LR

A["v₁"]
B["v₂"]
C(("v₃"))
D["v₄"]

A & B --> C --> D

The functions corresponding to nodes $v_3$ and $v_4$ are as follows:
$f_{v_3}(x)=2x_1+x_2\text{and}{f}_{v_4}(x)=\text{relu}(x)$
The formulations for each of the corresponding nodes will be:

$\text{Node }{v}_3:F_{v_3}\triangleq v_3^o=2v_3^{in,1}+v_3^{in,1}$
$\text{Node }{v}_4:F_{v_4}\triangleq (v_4^{in,1}>0⟹v_4^o=v_4^{in,1})\wedge (v_4^{in,1}\le 0⟹v_4^o=0)$

$$\begin{array}{r}{F}_{o\to v_3}\triangleq (v_3^{in,1}=v_1^o)\wedge (v_3^{in,2}=v_2^o)\\ F_{o\to v_4}\triangleq (v_4^{in,1}=v_3^o)\end{array}$$

Encoding of $G$ : $F_G\triangleq \underset{F_V}{\underbrace{F_{v_3}\wedge F_{v_4}}}\wedge \underset{F_E}{\underbrace{F_{o\to v_3}\wedge F_{o\to v_4}}}$
So far we have been considering the linear operations that form the weights and biases for the nodes in a neural network, however, there are also non-linearities that are present.

How to deal with non-linear activations?

We can create representations of non-linear activation functions by over-approximation of said functions. Over-approximation is the partitioning of the function into different input domains and binding the outputs corresponding to each of these input domains to some defined output domain. For example; we can make an over-approximation for the Sigmoid ($\sigma (x)$) function by diving the function into appropriate intervals that roughly corresponds to its overall behaviour, i.e., anything less than a certain threshold will give us a 0, anything more than a certain threshold will give us 1 and any input values in between them will give us a value of roughly 0.5.

$$\begin{array}{r}{F}_v\triangleq (v^{in,1}\le -1⟹0\le v^o\le 0.26)\wedge \\ (-1<v^{in,1}\le 0⟹0.26<v^o\le 0.5)\wedge \\ (0<v^{in,1}\le 1⟹0.5<v^o\le 0.73)\wedge \\ (v^{in,1}>1⟹0.73<v^o\le 1)\end{array}$$

The above can be generalized to any monotonically increasing or decreasing function $f_v$ (which is all activation functions). Assume $f_v$ is monotonically increasing, sample a sequence of real values $c_1<\cdots <c_n$

$$\begin{array}{r}{F}_v\triangleq (v^{in,1}\le c_1⟹lb<v^o\le f(c_2))\\ \wedge (c_1<v^{in,1}\le c_2⟹f_v(c_1)<v^o\le f(c_2))\\ ⋮\\ \wedge (c_n<v^{in,1}⟹f_v(c_n)<v^o\le ub\end{array}$$

where $lb$ and $ub$ are the *lower bound and upper bound respectively. A point to note in the case of over-approximations are that they will give us soundness, but not completeness. Completeness essentially means that our encoding can find counterexamples, which we are abandoning in this case. Soundness means being able to prove correctness properties using an encoding, which is of highest priority for our models.

A Concrete Example for checking Robustness

$$f_G:\mathbb{R}\to {\mathbb{R}}^2$$

Where the output vectors in ${\mathbb{R}}^2$ correspond to:

$$\left[\begin{array}{c}\\ p(\text{cat})\\ \\ p(\text{dog})\\ \end{array}\right]$$

Let the specification for this be:

$$\begin{array}{r}\{|x-c|\le 0.1\}\\ r\leftarrow f_G(x)\\ \{r_1>r_2\}\end{array}$$

where $c$ is the original image and $x$ is the perturbed image. The formula generated to check this statement is called a Verification Condition (VC), which if valid then ensures that the correctness property holds.

$$(\underset{p}{\underbrace{\text{precondtion}}}\wedge \underset{G}{\underbrace{\text{neural network}}})⟹\underset{q}{\underbrace{\text{postcondition}}}$$

inputs: $\{v_1,\dots ,v_n\}$ outputs: $\{v_{n+1},v_{n+2}\}$ Neural Network: $F_G$

$$\underset{\text{precondition}}{\underbrace{\left(\underset{i=1}{\overset{n}{\bigwedge }}|x_i-c_i|\le 0.1\right)}}\wedge \underset{\text{network}}{\underbrace{F_G}}\wedge \underset{\text{network input}}{\underbrace{\left(\underset{i=1}{\overset{n}{\bigwedge }}{x}_i=v_i^o\right)}}\wedge \underset{\text{network output}}{\underbrace{(r_1=v_{n+1}^o\wedge r_2=v_{n+2}^o)}}⟹\underset{\text{postcondition}}{\underbrace{(r_1>r_2)}}$$

LRA does not support vector operations hence the vector operations are decomposed into their constituent scalars. Absolute value operators are also not available in LRA, so we encode them as: $|x|\le 5\stackrel{LRA}{\to }(x\le 5)\wedge (-x\le 5)$. We need to connect variables of $F_G$ with inputs $x$ and output $r$.

For multiple Neural Networks, Correctness Properties have the form:

$$\begin{array}{r}\{P\}\\ r_1\leftarrow f_{G_1}(x_1)\\ r_2\leftarrow f_{G_2}(x_2)\\ ⋮\\ r_l\leftarrow f_{G_l}(x_l)\\ \{Q\}\end{array}$$

For the above specification we have:

$$\left(P\wedge \underset{i=1}{\overset{l}{\bigwedge }}{F}_i\right)⟹Q$$

For the above $F_i$ corresponds to $r_i\leftarrow f_{G_i}(x_i)$. $F_i$ combines encoding of the neural network $F_{G_i}$ along with connections with inputs and outputs, $x_i$ and $r_i$, respectively.

$$F_i\triangleq F_{G_i}\wedge \left(\underset{j=1}{\overset{n}{\bigwedge }}{x}_{i,j}=v_i^o\right)\wedge \left(\underset{j=1}{\overset{m}{\bigwedge }}{r}_{i,j}=v_{n+j}^o\right)$$ $$x_{i,j}\to v_i^o\to v_i^{in,1}$$

(soundness) If $F$ is valid, then the correctness property is true. (completeness) If it is invalid, we know there is a model $I\models \neg F$ if $F$ is encodeable in LRA.

Assuming, input and output variables of the encoding of $G_i$ are $v_1,\dots ,v_n$ and $v_{n+1},\dots ,v_{n+m}$; each graph $G_i$ has unique nodes and therefore input/output variables.

$$\left(P{\wedge }_{i=l}^l{F}_i\right)⟹Q$$

The above reads as: “If the precondition is true and we execute all $l$ networks, then the postcondition should be true”

Example:

$$\begin{array}{r}\{|x-1|\le 0.1\}\\ r\leftarrow f(x)\\ \{r\ge 1\}\end{array}$$

where $f(x)=x$.

We take the input value $x$ to be equal to 0.99: $x=0.99$. This is a valid value of $x$ as per the defined precondition.
Therefore, we have: $f(x)⟹f(0.99)=0.99$
$F$ is invalid as $I\models \neg F$, thus the correctness property does not hold.

Info

In the formulations that we have gone over, disjunctions arise due to the ReLU function, due to its active/inactive states for all possible inputs to the Neural Network. Disjunctions cause problem during solving, as without them LRA (and a similar system called MILP) are polynomial time solvable. To solve LRAs with disjunction we either simplify the formulations by using lightweight techniques to discover whether ReLUs are active or not (abstraction-based verification). Alternatively, we can add additional bias to make all the ReLUs either always active or always inactive. Another thing to consider is that verified NNs in LRA may not really be robust when considering bit-level behaviour.

DPLL (Davis-Putnam-Logemann-Loveland)

DPLL algorithm checks the satisfiability of Boolean formulae and underlies modern SMT and SAT solver. An extension of DPLL algorithm is needed to handle first-order formulae over different theories.

Conjunctive Normal Form (CNF)

DPLL expects formulae that will be inputted to be in the shape of Conjunctive Normal Form (CNF), all Boolean formulae must be written in CNF as can be seen below

$$\text{CNF: }{C}_1\wedge \cdots \wedge C_n$$

where each sub-formula $C_i$ is called a clause and is for the form

$$l_1\vee \cdots \vee l_{m_i}$$

where each $l_i$ is called a literal and is either a variable ($p$) or its negation ($\neg p$). Thus the structure of an entire input should be in the form:

flowchart TD

C_1((C₁))
C_2((C₂))
C_n((Cₙ))
l_1((l₁))
l_2((l₂))
l_m((lₘ))
p((p))
np((¬p))
d1([......])
d2([......])

C_1 --> l_1 & l_2 & l_m
l_1 --> p & np


subgraph Clauses
C_1 -. ∧ .- C_2 -. ∧ .- d1 -. ∧ .- C_n
end


subgraph Literals
l_1 -. ∨ .- l_2 -. ∨ .- d2 -. ∨ .- l_m
end

subgraph Variables
p
np
end

DPLL has the following two alternating phases

flowchart LR

A(["Deduction"])
B(["Search"])

A --> B --> A

Deduction

Boolean Constant Propagation (BCP)

The algorithm searches the Boolean CNF for clauses with single literals, i.e., clauses that contain only one Boolean variable instead of disjunctions of multiple literals. For a model to be satisfiable, the single literal must be $\text{True}$ if the variable is $p$ or $\text{False}$ if the variable is $\neg p$.

$$(l)\wedge C_2\wedge \cdots \wedge C_n$$

The BCP phase will look for all unit clauses and replace their literals with $\text{True}$.

Example:

$$F\triangleq (p)\wedge (\neg p\vee r)\wedge (\neg r\vee q)$$ $$\begin{array}{rl}\text{BCP: }& (\text{True})\wedge (\neg \text{True}\vee r)\wedge (\neg r\vee q)\\ & \equiv (r)\wedge (\neg r\vee q)\\ & \\ \text{BCP: }& (\text{True})\wedge (\neg \text{True}\vee q)\\ & \equiv q\\ & \\ \text{BCP: }& (\text{True})\end{array}$$

Thus $F$ is SAT with the model

$$\{p\mapsto \text{True},q\mapsto \text{True},r\mapsto \text{True}\}$$

Deduction + Search

DPLL:
	Data: A formula F in CNF
	Result: I ⊧ F or UNSAT
▹ Boolean Constant Propagation (BCP)
	while there is a unit clause (l) in F do:
		Let F be F[F ↦ True]
	if F is True then return SAT
▹ Search
	for every possible variable in F do:
		If DPLL(F[p ↦ True]) is SAT then return SAT
		If DPLL(F[p ↦ False]) is SAT then return SAT
	return UNSAT

The model $I$ that is returned by DPLL when the input is SAT is maintained implicitly in the sequence of assignments to variables (of the form $[l\mapsto \dots ]$ and $[p\mapsto \dots ]$).

$$\begin{array}{rl}F\triangleq & (p\vee r)\wedge (\neg p\vee q)\wedge (\neg p\vee \neg r)\\ \text{recursion 1: }& F_1=F[p\mapsto \text{True}]=q\wedge (\neg q\vee \neg r)\\ \text{recursion 2: }& F_2=F_1[q\mapsto \text{True}]=(\neg r)\\ & F_3=F_2[p\mapsto \text{FAlse}]=(\text{True})\end{array}$$

DPLL returns SAT and then implicitly builds the model for $F$:

$$\{p\mapsto \text{True},q\mapsto \text{True},r\mapsto \text{False}\}$$

Partial Model: DPLL can terminate with SAT and without assigning values to each and every variable, these incomplete models are called as partial models. The unfilled variables are essentially don’t care variables and can be filled in any way we want. For $F\triangleq p\wedge (q\vee p\vee \neg r)\wedge (p\vee \neg q)$, $I=\{p\mapsto \text{True}\}$ is a partial model.

DPLL Modulo Theories

DPLL modulo theories or DPLLᵀ are extensions of DPLL over formulae in mathematical theories such as LRA. We essentially treat a formula completely by taking it as a Boolean, then incrementally add more and more theory info to conclusively prove SAT or UNSAT.

$$\begin{array}{rl}F\triangleq & (x\le 0\vee x\le 10)\wedge (\neg x\le 0)\\ & \\ F^B\triangleq & (p\vee q)\wedge (\neg p)\end{array}$$

flowchart LR

A[Original Formula
F]
B[Boolean Abstraction
Fᴮ]

A --B---> B --T---> A

DPLLᵀ Algorithm

Constraints are lost in the abstraction process ($F\stackrel{B}{\to }{F}^B$) such as the relation between the different inequalities. If $F^B$ is UNSAT, then $F$ is UNSAT, but, if $F^B$ is SAT, it does not imply that $F$ is SAT.
The DPLLᵀ algorithm first uses DPLL to check if $F^B$ is UNSAT following the properties of abstraction. IF $F^B$ is SAT, we will take the model $I$ returned by DPLL($F^B$) and map it to the formula $I^T$ in the theory we have abstracted the theory from (LRA in our case). If the theory solver deems $I^T$ satisfiable, $F$ is satisfiable. Otherwise, DPLLᵀ learns that $I^T$ is not a model so it negates $I$ and conjoins it to $F^B$. DPLLᵀ lazily learns more and more facts about the formula and refines the abstraction until the algorithm can decide SAT or UNSAT. The whole process is illustrated below:

flowchart LR

FU([F is UNSAT])
FB((Fᴮ))
I((I))
IT([Iᵀ])
FS([F is SAT])
NI((¬I))

FB --SAT---> I --T--> IT --SAT---> FS
FB --UNSAT---> FU
IT --UNSAT---> NI --∧---> FB

DPLL takes care of the disjunctions by taking clauses and searching for satisfiability by mapping literals to $\text{True}$ or $\text{False}$ values. DPLL assumes access to a theory solver to take care of the conjunctions of linear inequalities, in order to check their satisfiability. For LRA, the Simplex Algorithm can be used as a theory solver.

Example:

$$\begin{array}{rl}\text{LRA: }& (x\ge 10)\wedge ((x<0)\vee (y\ge 0))\\ F^B\text{: }& p\wedge (q\vee r)\\ & \\ I_1=& \{p\mapsto \text{True},q\mapsto \text{True}\}\\ & p\wedge q\text{ : SAT}\\ & \\ I_1^T=& \underset{p}{\underbrace{x\ge 10}}\wedge \underset{q}{\underbrace{x<0}}\text{ : UNSAT}\\ & \\ F^B\wedge \neg I_1:& p\wedge (q\vee r)\wedge \underset{\neg I_1}{\underbrace{(\neg p\vee \neg q)}}\\ I_2=& p\wedge \neg q\wedge r\\ & \\ I_2^T:& \{x\mapsto 10,y\mapsto 0\}\text{ : SAT}\end{array}$$

How to convert Boolean formulae into CNF? Usually we use De Morgan’s Law (however all are $O(exp)$) What can be used: Tseitin’s Transformation ($O(n)$)

Tseitin’s Transformation

$$\underset{x\in S_{\text{var}}}{F}\stackrel{\text{Tseitin’s Transformation}}{\to }\underset{x\in S_{\text{var}}^{\prime }}{F^{\prime }}\text{(CNF)}$$ $$S_{\text{var}}\subseteq S_{\text{var}}^{\prime }$$

Any model of $F^{\prime }$ is also a model of $F$, if we disregard the interpretations of newly added variables. If $F^{\prime }$ is UNSAT, then $F$ is UNSAT, so we just need to invoke DPLL on $F^{\prime }$.

Tseitin’s transformation changes a formula of computations into a set of instructions, each containing one or two variables, connected by a single unary or binary operator respectively. For example:

def f(x, y, z):
    pass
    return x + (2*y + 3)

The above function can be decomposed into a set of instructions, that when executed sequentially will provide the same result. For the above example, the function can be decomposed into the following instructions:

def f(x,y,z):
    t1 = 2 * y
    t2 = t1 + 3
    t3 = x + t2
    return t3

If we are able to conjoin these sequential set of instructions, we will be converting the complex computation in non-CNF to CNF with added temporary variables $t_i$; Tseitin’s transformation follows a similar procedure over Boolean formulae to generate a CNF.

Tseitin Step 1: NNF

Negation Normal Form (NNF) is achieved by pushing negation inwards so that $\neg$ only appears next to variables, e.g., $\neg p\vee \neg r$ instead of $\neg (p\wedge r)$.

$$\begin{array}{rl}\neg (F_1\wedge F_2)& =\neg F_1\vee \neg F_2\\ & \\ \neg (F_1\vee F_2)& =\neg F_1\wedge \neg F_2\\ & \\ \neg \neg F_1& =F_1\end{array}$$

Tseitin Step 2: Subformula Rewriting

Any subformula of $F$ that contains a conjunction/disjunction is called a subformula. We don’t consider subformula at literal level.

$$F\triangleq \underset{F_4}{\underbrace{(\underset{F_1}{\underbrace{p\wedge q}})\vee (\underset{F_3}{\underbrace{\underset{F_2}{\underbrace{q\wedge \neg r}}\wedge s}})}}$$

$F_1,F_2$ are the deepest level of nesting, $F_2$ is subformula of $F_3$. And all $F_i$ are subformulae of $F_4$.

Assuming $F$ has $n$ subformulae:

• For every subformula $F_i$ of $F$, create a fresh variable $t_i$. These variables are analogous to the temporary variables that was introduced to decompose some complex program
• Starting with the most deeply-rooted subformula: let $F_i$ be of the form $l_i\circ l_i^{\prime }$, where $\circ$ is $\wedge$ or $\vee$ and $l_i,l_i^{\prime }$ are literals. One or both of $l_i$ and $l_i^{\prime }$ may be the new variable $t_j$ denoting a subformula $F_j$ of $F_i$, create the formula:

$$F_i^{\prime }\triangleq t_i\iff (l_i\circ l_i^{\prime })$$

These formulae are analogous to the assignments to temporary variables in code, where $\iff$ is the logical analogue of variable assignment (=).

$$\begin{array}{rl}{F}_1^{\prime }& \triangleq t_1\iff (p\wedge q)\\ F_2^{\prime }& \triangleq t_2\iff (q\wedge \neg r)\\ F_3^{\prime }& \triangleq t_3\iff (t_2\wedge s)\\ F_4^{\prime }& \triangleq t_4\iff (t_1\vee t_3)\end{array}$$ $$\begin{array}{rl}{l}_1\iff (l_2\vee l_3)& \equiv (\neg l_1\vee l_2\vee l_3)\wedge (l_1\vee \neg l_2)\wedge (l_1\vee \neg l_3)\\ l_1\iff (l_2\wedge l_3)& \equiv (\neg l_1\vee l_2)\wedge (\neg l_1\vee l_3)\wedge (l_1\vee \neg l_2\vee \neg l_3)\end{array}$$ $$F^{\prime }\triangleq t_n\underset{i}{\bigwedge }{F}_i^{\prime }$$

each $t_i$ is assigned true iff subformula $F_i$ evaluates to true. The constant $t_n$ in $F^{\prime }$ says that $F$ must be true, $t_n$ is like a return statement.

$$F^{\prime }\triangleq t_n\wedge F_1^{\prime }\wedge F_2^{\prime }\wedge F_3^{\prime }\wedge F_4^{\prime }$$

Theory Solving

Theory solver for LRA receives $F$ formula as a conjunction of linear inequalities

$$\underset{i=1}{\overset{n}{\bigwedge }}\left(\sum _{j=1}^m{c}_{ij}\cdot x_j\ge b_i\right)\text{ where }{c}_{ij},b_i\in \mathbb{R}$$

Goal: Check SAT for $F$ and discover $I\models F$.

The Simplex Algorithm will be used as a theory solver, and it expects formulae to be conjunctions of equalities of the form ${\sum }_i{c}_i\cdot x_i=o$ and bounds of the form $l_i\le x_i\le u_i$, $u_i,l_i\in \mathbb{R}\cup \{\infty ,-\infty \}$. The infinities are included to ensure that variables with no upper/lower bounds are also adequately represented.

Converting inequalities into simplex (slack) form:

$$\begin{array}{c}F\triangleq \underset{i=1}{\overset{n}{\bigwedge }}\left(\sum _{j=1}^m{c}_{ij}\cdot x_j\ge b_i\right)\\ s_i=\sum _{j=1}^m{c}_{ij}\cdot x_j\text{ AND }{s}_i\ge b_i\end{array}$$

where $s_i$: Slack Variable (analogous to tseitin temporary variables)

Let $F_s$ be the simplex form of some formula $F$. Then we have the following guarantees (analogue of Tseitin Transformation for non-CNF formulae):

• Any model of $F_s$ is a model of $F$, disregarding assignments to slack variables.
• If $F_s$ is UNSAT, then $F$ is UNSAT.

Simplex Algorithm

Goal: Find a satisfying argument that maximizes some objective function. Our interest in verification is to find any satisfying assignment, so it will be a subset of Simplex.

The set of variables in simplex form is classified into two subsets:

• Basic Variables: those that appear on the left hand side of the equality; initially basic variables are the slack variables.
• Non-Basic Variables: all other variables.

At the beginning, basic variables: $\{s_1,s_2,s_3\}$ and non-basic: $\{x,y\}$, as Simplex progresses, formulae are rewritten so some basic variables may become non-basic and vice versa.

Simplex simultaneously looks for a model and a proof for unsatisfiability. Each equality defines a halfspace (something that splits ${\mathbb{R}}^2$ into two parts). Simplex starts from $I_0=\{x\mapsto 0,y\mapsto 0\}$ and toggles the values to satisfy all of the equalities. We are total ordering our variables to make it easier to refer to specific variables. We assumes variables are of the form $x_1,\dots ,x_n$. Given a basic variable $x_i$ and a non-basic variables $x_j$, we will use $c_{ij}$ to denote the coefficient of $x_j$ in the equality

$$x_j=\cdots +c_{ij}\cdot x_j+\dots$$

for variable $x_i$, $l_i$ and $u_i$ denotes its lower and upper bound. (Non-Slack variables have no bounds)

Two invariants that are maintained in Simplex:

• $I$ always satisfies the equalities, so only bounds may be violated; initially true as all $0$ in $I$.
• Bounds of all non-basic variables are satisfied; initially true as they have no bounds.

For $x_i<l_i$, we need to update $x_i$ using any non-basic var $x_j$. Pick any $x_j$ with $c_{ij}\ne 0$. If no such $x_j$ exists that satisfies our condition $⟹$ UNSAT Increase current interpretation by $\frac{l_i-I(x_i)}{c_{ij}}$, interpretation of $x_j$ increases by $l_i-I(x_i)$, barely satisfying the lower bound, i.e., $I(x_i)=l_i$.

$$\begin{array}{rlrl}& x_j\mapsto 0& & x_j\mapsto \frac{l_i-I(x_i)}{c_{ij}}\\ & x_i=\cdots +0\dots & & x_i=\cdots +c_{ij}\cdot \frac{l_i-I(x_i)}{c_{ij}}\dots \\ & & & x_i=\cdots +l_i-I(x_i)\dots \\ & \because I(x_i)=0\text{ for previous interpretation}& & \\ & & & x_i=\cdots +l_i+\dots \\ & & I(x_i)=\underset{\text{previous}}{\underbrace{I(x_i)}}+l_i& \end{array}$$

After updating $x_j$, there is a chance that we might have violated the bounds of $x_j$ so we rewrite the equation such that $x_j$ becomes the basic variable and $x_i$ becomes the non-basic variable. The pivot operation follows as:

$$\begin{array}{c}{x}_i=\sum _{k\in N}{c}_{ik}{x}_k\\ \\ x_j=\underset{\text{replace }{x}_j\text{ with this}}{\underbrace{-\frac{x_i}{c_{ij}}+\sum _{k\in N/\{j\}}\frac{c_{ik}}{c_{ij}}{x}_k}}\end{array}$$

where $N$ is the set of indices of non-basic variables. For Simplex basic variables are dependent variables and non-basic variables are independent variables.

Example:

$$\begin{array}{c}x+y\ge 0\\ -2x+y\ge 2\\ -10x+y\ge -5\end{array}$$

Converting the set of equations into their slack form:

$$\begin{array}{c}{s}_1=x+y\\ s_2=-2x+y\\ s_3=-10x+y\\ s_1\ge 0\\ s_2\ge 2\\ s_3\ge -5\end{array}$$

Ordering $x,y,s_1,s_2,s_3$ and apply first model:

$$I_0=\{x\mapsto 0,y\mapsto 0,s_1\mapsto 0,s_2\mapsto 0,s_3\mapsto 0\}$$

$s_1$ and $s_3$ are satisfied, $s_2$ is not: $I_0(s_2)=0$ but $s_2\ge 2$.

Modulate the first element in our ordering: $x$ Decrease $I_0(x)$ to $-1$ to meet $s_2$‘s bounds.

$$I_1=\{x\mapsto -1,y\mapsto 0,s_1\mapsto -1,s_2\mapsto 2,s_3\mapsto 10\}$$

Now we pivot with unchanged bounds:

$$\begin{array}{c}x=0.5y-0.5s_2\\ s_1=1.5y-0.5s_2\\ s_3=-4y+5s_2\end{array}$$

$I_1(s_1)=-1<0$, bound not satisfied.

New ordering of variables: $\underbrace{y},s_2,x,s_1,s_3$ $I(y)$ value increased by $\frac{1}{1.5}=2.3$

$$I_2=\left\{x\mapsto -\frac{2}{3},y\mapsto \frac{2}{3},s_1\mapsto 0,s_2\mapsto 2,s_3\mapsto \frac{7}{3}\right\}$$

At this point we pivot $y$ with $s_1$

Simplex terminates as $I_2\models F$.

Simplex terminates due to the fact that the variables are ordered and we always look for the first variables violating bounds (Bland’s Rule), this ensures that we never revisit the same set of basic and non-basic variables.

Note

Basic variables are dependent variables and Non-Basic variables are independent variables.

Using Simplex as the theory solver within DPLLᵀ allows us to solve for LRA, but this approach is not scalable as ReLUs are encoded as disjunctions. This is because the SAt-solving part of DPLLᵀ will handle it and consider every possible case of disjunction (active = x, inactive = 0), leading to many calls to Simplex, so we extend Simplex to Reluplex.

Reluplex

Equations in reluplex form:

1. equations (same as simplex)
2. bounds (same as simplex)
3. ReLU constraints of the form: $x_i=relu(x_j)$ (also add bound(s) $x_i\ge 0$ implied by relu)

We call simplex on the weaker version (less constrained) formula of $F$ called $F^{\prime }$, which does not have any relu constraints. If $F^{\prime }$ is UNSAT, then $F$ is UNSAT.

$$F⟹F^{\prime }\text{ (is valid)}$$

If Simplex returns $I\models F^{\prime }$, it may not be that $I$ satisfies $F$. If $I/\models F$, we pick one of the violated relu constraints

$$x_i=relu(x_j)$$

and modify $I$ to make sure that it is not violated. If any of $x_i$ or $x_j$ are basic, we pivot it with a non-basic variable. The pseudocode for Reluplex Algorithm can be found below:

Note: Without Case Splitting, reluplex might not terminate - it may get stuck in a loop where the Simplex satisfies all bounds but violates a relu, then satisfying that relu causes a bound to be violated and so on.

We try to count whether a relu constraint has not been attempted to be fixed for more than $\tau$ times. $x_i=relu(x_j)$ is split into two cases: $F_1\triangleq x_j\ge 0\wedge x_i=x_j$ and $F_2\triangleq x_j\le 0\wedge x_i=0$. Reluplex is recursively invoked on the two instances of the problem: $F\wedge F_1$ and $F\wedge F_2$.

If any of the instances are SAT, then $F$ is SAT.

$$F\equiv (F\wedge F_1)\vee (F\wedge F_2)$$

Abstraction Based Verification

Approximate (or Abstract) techniques for verification can have two possible outcomes:

• If they succeed, they can produce proofs of correctness for a neural network.
• If they fail, we do not know whether the correctness property holds or not.

Abstraction Based Verification techniques require specific techniques in order to abstractly define a domain of inputs to be given to a modified version of a neural network in order to test their robustness. There are mainly three ways of abstracting input domains for verification: Interval Abstraction, Zonotope Domain Abstraction and Polyhedra Domain Abstraction.

flowchart TD

A["Abstraction Techniques"] --> B["Interval"] & C["Zonotope"] & D["Polyhedra"]

In Interval Abstraction, we define an interval for each of the entire set of variables/predicates independently of each other, thereby defining a rectangular structure in some high dimensional space (hyperrectangles).

In Zonotope Domains, instead of defining hyperrectangles, the input domain is defined as a collection of parallelograms superimposed on each other created using auxiliary variables called generators whose values range from -1 to 1, each generator has its corresponding coefficient values which help in making a convex structure. The main difference between zonotopes and interval domains is that we can encode the relations between multiple variables/predicates using the generators.

In Polyhedra Domain, we have a setup similar to that of Zonotope Domains, however instead of the generator values ranging from -1 to 1, we have a set of linear inequalities bound together via conjunctions that determine the overall convex shape of the interval we are defining.

These techniques have been described in more detail below:

Neural Interval Abstraction

Using the specifications that have been defined beforehand, we see that it defines a system to test for individual images through a particular neural network.

$$\begin{array}{r}\{|x-c|\le 0.1\}\\ r\leftarrow f(x)\\ \{class(r)=1\}\end{array}$$

There is a great number of possible images $x$ that can be tested for, so we will try to lift the function $f(\dots )$ to work over a set of images instead in order to comprehensively test for possible examples. Thus we want to transform the neural network computation that is represented by $f$

$$f:{\mathbb{R}}^n\to {\mathbb{R}}^m$$

to

$$f^S:\mathcal{P}({\mathbb{R}}^n)\to \mathcal{P}({\mathbb{R}}^m)$$

where $\mathcal{P}(S)$ is the power set of set $S$. Thus, the newly defined mapping from a set of elements ($X$) to the set of all predictions corresponding to all elements of $X$ will following form:

$$f^S(X)=\{y|x\in X,y=f(x)\}$$

where $X$ is a set that is defined as per the specifications as the set of all $x\in X$ that are valid under the precondition. Which for our example will be:

$$X=\{x||x-c|\le 0.1\}$$

To verify our property, we simply check:

$$f^S(X)\subseteq \{y|class(y)=1\}$$

all runs of $f$ on every image $x\in X$ result in network predicting class $1$. We are defining $\infty$ set of inputs using data structures that we can manipulate called abstract domains. So we are taking a Neural Network and generating a version that can take a potentially infinite set of images, however during this abstraction process we will also be losing some precision, as we will be seeing soon. $f^S(X)$ is called the concrete transformer of $f$.

Example: The function is: $f(x)=x+1$ Corresponding concrete transformer: $f^S(X)=\{x+1|x\in X\}$

Interval Abstract Domain considers an interval over $\mathbb{R}:[l,u]$ where $l,u\in \mathbb{R}$ and $l\le u$

$$\{x|l\le x\le u\}$$

We simplify concrete transformer by considering only sets which have a nice form (Abstract Interpretation).

$$f^a\left(\underset{\text{input interval}}{\underbrace{[l,u]}}\right)=\underset{\text{output interval}}{\underbrace{[l+1,u+1]}}$$

$f^a$ is the abstract transformer of $f$. Interval $[l,u]$ is infinite (considering $l<u$) so $f^a$ adds $1$ to an infinite set of reals. Expanding the above to any arbitrary $n^{th}$ dimension we have a n-dimensional interval (hyperrectangle) region in ${\mathbb{R}}^n$, i.e., set of all n-ary vectors $\{\overrightarrow{x}\in {\mathbb{R}}^n|l_i\le x_i\le u_i\}$

$$\begin{array}{rlr}& {\mathbb{R}}^n\text{ : hyperrectangle interval : }& ([l_1,u_1],\dots ,[l_i,u_i],\dots ,[l_n,u_n])\\ & \uparrow & \\ & ⋮& \\ & {\mathbb{R}}^3\text{ : box interval : }& ([l_1,u_1],[l_2,u_2],[l_3,u_3])\\ & \uparrow & \\ & {\mathbb{R}}^2\text{ : rectangle interval : }& ([l_1,u_1],[l_2,u_2])\\ & \uparrow & \\ & \mathbb{R}\text{ : line interval : }& ([l,u])\end{array}$$

Soundness: We need to ensure that the $f^a$ designed is a sound approximation of $f^S$. Output of $f^a$ should be a superset of $f^S$ , as we do not want to miss any behaviour. For any interval $[l,u]$, we have: $f^S([l,u])\subseteq f^a([l,u])$. Equivalently, for any $x\in [l,u]$, we have: $f(x)\in f^a([l,u])$ Although in practice we often see: $f^S([l,u])\subset f^a([l,u])$

We are modifying the functions to take in intervals of inputs and output the interval that contains the set of all the mappings from the input interval. However, the interval domain cannot capture the relations between different dimensions (non-relational)

$$X=\{(x,x)|0\le x\le 1\}$$

Best we can do is define the unit square between $(0,0)$ and $(1,1)$ which is denoted as the 2D interval $([0,1],[0,1])$. $X$ defines points where higher x-coordinates correspond to higher y-coordinates. Our abstract abstract domain can only represent rectangles whose faces are parallel to the axes. So instead of capturing the relation between two dimensions, we are simply saying that any value that is in $[0,1]$ for $x$ can correspond to any value of $y$ in $[0,1]$, thereby overapproximating for the function to the extent that we are unable to conserve any information about the relations between the two elements of the input.

Example 1: Consider $f(x,y)=x+y$

$f^S:\mathcal{P}({\mathbb{R}}^2)\to \mathcal{P}(\mathbb{R})$ is defined as $f^S(X)=\{x+y|(x,y)\in X\}$
$f^a$ is defined as a function that takes two intervals (a rectangle) representing range of values for $x_1$ and $x_2$

$$\begin{array}{c}{f}^a([l,u],[l^{\prime },u^{\prime }])=[l+l^{\prime },u+u^{\prime }]\\ \\ f^a([1,5],[100,200])=[101,205]\end{array}$$

Take any $(x,y)\in ([l,u],[l^{\prime },u^{\prime }])$ By definition, $l\le x\le u$ and $l^{\prime }\le y\le u^{\prime }$, so $l+l^{\prime }\le x+y\le u+u^{\prime }$ thus $x+y\in [l+l^{\prime },u+u^{\prime }]$ proving soundness.

Example 2: Consider $f(x,y)=x\times y$ For only positive inputs to $f^a$ : $f^a([l,u],[l^{\prime },u^{\prime }])=[l\times l^{\prime },u\times u^{\prime }]$ However,

$$f^a([-1,1],[-3,-2])=[3,-2]$$

We were expecting $l\times l^{\prime }\le u\times u^{\prime }$ but we got $l\times l^{\prime }\ge u\times u^{\prime }$ which would mean that the result is not a proper interval. Thus we have to compute the interval domain in a different way:

$$f^a([l,u],[l^{\prime },u^{\prime }])=[min(B),max(B)]$$

where

$$B=\{l\times l^{\prime },l\times u^{\prime },u\times l^{\prime },u\times u^{\prime }\}$$

For the above numerical example we therefore have: $f^a([-1,1],[-3,-2])=[min(B),max(B)]=[-3,3]$

These were all basic abstract transformers that overapproximate for some given binary function.

Affine Function:

$f(x_1,\dots ,x_n)={\sum }_i{c}_i{x}_i$ where $c_i\in \mathbb{R}$, we can define the abstract transformer as:

$$f^a([l_1,u_1],\dots ,[l_n,u_n])=\left[\sum _i{l}_i^{\prime },\sum _i{u}_i^{\prime }\right]$$

where $l_i^{\prime }=\text{min}(c_i{l}_i,c_i{u}_i)$ and $u_i^{\prime }=\text{max}(c_i{l}_i,c_i{u}_i)$, trying to cover the largest area possible.

Example: $f(x,y)=3x+2y$ $f([5,10],[20,30])=[3\times 5+2\times 20,3\times 10+2\times 30]=[55,90]$

Monotonic Function:

$f:\mathbb{R}\to \mathbb{R}$ $f^a([l,u])=[f(l),f(u)]$

Example: ReLU in the domain $[3,5]$ $relu(3)\le relu(x)\le relu(5)$ Therefore, $f^a([3,5])=[relu(3),relu(5)]$

Composing Abstract Transformers:

For $(f\circ g)(x)$ we do not need to define a different transformer, we can define one for $f(f^a)$ and one for $g(g^a)$ and compose them to find a sound abstract transformer of $f\circ g$ : $f^a\circ g^a$

Example: $g(x)=3x$, $f(x)=relu(x)$ and $h(x)=f(g(x))$ Affine followed by ReLU and output in $\mathbb{R}$

$$h^a([2,3])=f^a(g^a([2,3]))=f^a([6,9])=[6,9]$$

Abstractly Interpreting Neural Networks

For a network $G=(V,E)$ we have: $f_G:{\mathbb{R}}^n\to {\mathbb{R}}^m$ where $n=|V^{in}|,m=|V^{out}|$

• $f_G^a$ takes $n$ intervals and outputs $m$ intervals
• for every output node $v_i$: ${\text{out}}^a(v_i)=[l_i,u_i]$ (note: we have a fixed ordering of nodes)
• for every non-input node $v$ : ${\text{out}}^a(v)=f_v^a({\text{out}}^a(v_1),\dots ,{\text{out}}^a(v_k))$ where $f_v^a$ is the abstract transformer of $f_v$ and $v$ has incoming edges $(v_1,v),\dots ,(v_k,v)$
• the output of $f_v^a$ is the set of intervals ${\text{out}}^a(v_1),\dots ,{\text{out}}^a(v_m)$ where $v_1,\dots ,v_m$ are output nodes

Example: $f_{v_3}(x)=2x_1+x_2$ $f_{v_4}(x)=relu(x)$

$f_G^a([0,1],[2,3])$:

flowchart LR

v1((v₁))
v2((v₂))
v3((v₃))
v4((v₄))
o["[2,5]"]
v1 --[0,1]--> v3
v2 --[2,3]--> v3
v3 --[2,5]--> v4 --> o

$$\begin{array}{c}{\text{out}}^a(v_1)=[0,1]\\ {\text{out}}^a(v_2)=[2,3]\\ {\text{out}}^a(v_3)=[2\times 0+2,2\times 1+3]=[2,5]\\ {\text{out}}^a(v_4)=[relu(2),relu(5)]\end{array}$$

Limitations

Interval Domain often overshoots and computes wildly overapproximated solutions.

Example 1:

flowchart LR

v1((v₁))
v2((v₂))
v3((v₃))
o["[-1,1]"]

v1 --[0,1]--> v2 --[-1,0]--> v3 --> o
v1 --[0,1]--> v3

$f_{v_3}(x)=-x$ $f_{v_4}(x)=x_1+x_2$

So we can see that $f_G(x)=0$ Expected abstract transformer should be: $f_G^a([l,u])=[0,0]\forall l,u\in \mathbb{R}$ But $f_G^a([l,u])$ returns $[-1,1]$ since it does not know the relation between $-x$ and $x$.

Example 2:

flowchart LR

v1((v₁))
v2((v₂))
v3((v₃))
o1["[0,1]"]
o2["[0,1]"]

v1 --[0,1]--> v2 & v3
v2 --> o1
v3 --> o2

$f_{v_2},f_{v_3}$ are both relus. $f_G(x)=(x,x)$ is the desirable result However, $f_G^a([0,1])=([0,1],[0,1])$ is the result. $f_G^a$ tells us that for inputs between $0$ and $1$, the neural network can output $(x,y)$ where $0\le x,y\le 1$ which is a loose approximation, we needed $0\le x\le 1$.

We cannot capture set of points where $x=y$ and $0\le x\le 1$, instead we can give $([0,1],[0,1])$, syntactically an abstract element in the interval domain is captured by constraints of the form:

$$\underset{i}{\bigwedge }{l}_i\le x_i\le u_i$$

Every inequality involves a single variable and fails to capture relationships so interval domain is called non-relational. To address these limitations we will be turning towards Zonotope-based abstraction methods.

Zonotope: Relational Abstract Domain

Assume that we have a set of $m$ real-valued generator variables ${ϵ}_1,\dots ,{ϵ}_m\in [-1,1]$. A $1D$ zonotope is the set of all points in the set

$$\left\{c_0+\sum _{i=1}^m{c}_i{ϵ}_i|{ϵ}_i\in [-1,1]\right\}$$

where $c_i\in \mathbb{R}$. For 1 generator variable:

$$\left\{c_0+c_1{ϵ}_1|{ϵ}_1\in [-1,1]\right\}$$

This is just defining an interval of the form $[c_0-c_1,c_0+c_1]$, assuming $c_1\ge 0$. In the defined interval, $c_0$ is the centre. For a one dimensional zonotope this can be interpreted as the point $c_0$ being stretched in both the available directions with the magnitude of $c_1$.

$$c_0-c_1\leftarrow c_0\to c_0+c_1$$

Zonotopes are more expressive from ${\mathbb{R}}^2$ and above. Thus, similarly for multiple generator variables, we can observe the same “stretching” effect. The initial coordinates are stretched into a line which represents a specific vector with endpoints at $(c_{10}+c_{11},c_{20}+c_{21})$ and $(c_{10}-c_{11},c_{20}-c_{21})$. The vector is then stretched by another generator into a parallelogram. The more generator variables there are, the more faces there are in the resulting zonotope. However each zonotope can be described a convex figure resulting from the summation of different parallelograms. The name zonotope is derived from zona meaning ‘belt’ in Latin. Zonotopes are named such as we can trace an uninterrupted path through the parallel vectors that wrap around the figure like a belt.

In n-dimensions, a zonotope with m-generators is the set of all points

$$\left\{\left(\underset{\text{first dimension}}{\underbrace{c_{10}+\sum _{i=1}^m{c}_{1i}{ϵ}_i}},\dots ,\underset{\text{nth dimension}}{\underbrace{c_{n0}+\sum _{i=1}^m{c}_{ni}{ϵ}_i}}\right)|{ϵ}_i\in [-1,1]\right\}$$

Example: $(1+{ϵ}_1,2+{ϵ}_2)$ $(1+{ϵ}_1+0{ϵ}_2,2+0{ϵ}_1+{ϵ}_2)$ Centre will be $(1,2)$, the centre of the zonotope is always the vector of the constant coefficients.

$(2+{ϵ}_1,2+{ϵ}_1)$ Two dimensions are equal, so we get a line shape centred at $(2,2)$. Zonotopes allow us to model relational intervals with the help of the generator variables. Above we see a zonotope modelling an interval that follows $(x,x)$.

$(2+{ϵ}_1,3+{ϵ}_1+{ϵ}_2)$ coefficients of ${ϵ}_1$ are $(1,1)$ so it stretches the centre $(2,3)$ along the $(1,1)$ vector

coefficients of ${ϵ}_2$ are $(0,1)$ so it stretches all points along $(0,1)$ vector

Above is the final zonotope, adding more and more faces adds more faces to the zonotope.

$$\left\{\left(\underset{\text{first dimension}}{\underbrace{c_{10}+\sum _{i=1}^m{c}_{1i}{ϵ}_i}},\dots ,\underset{\text{nth dimension}}{\underbrace{c_{n0}+\sum _{i=1}^m{c}_{ni}{ϵ}_i}}\right)|{ϵ}_i\in [-1,1]\right\}$$

We will use a compact notation to signify the above equation, it will be defined as a tuple of vectors of coefficients

$$(⟨c_{10},\dots ,c_{1m}⟩,\dots ,⟨c_{n0},\dots ,c_{nm}⟩)$$

for an even more compact notation

$$(⟨c_{1i}{⟩}_i,\dots ,⟨c_{ni}{⟩}_i)$$

where $i$ ranges from $0$ to $m$, the number of generators. We can compute the upper bound of the zonotope in the $j$th dimensionby solving the following optimization problem:

$$\begin{array}{c}\text{max }{c}_{j0}+\sum _{i=1}^m{c}_{ji}{ϵ}_i\\ \text{s.t. }{ϵ}_i\in [-1,1]\end{array}$$

which can easily be solved by setting ${ϵ}_i$ to $1$ if $c_{ji}>0$ and $-1$ otherwise.

Similarly, we can compute the lower bound of the zonotope in the j-th dimension by minimizing instead of maximizing, solving the optimization problem by setting ${ϵ}_i=-1$ if $c_{ji}>0$ and $0$ otherwise.

Example: $(2+{ϵ}_1,3+{ϵ}_1+{ϵ}_2)⟹(⟨2,1,0⟩,⟨3,1,1⟩)$

Upper Bound in vertical dimension: $3+{ϵ}_1+{ϵ}_2;{ϵ}_1,{ϵ}_2=1⟹3+1+1=5$

$f(x,y)=x+y$

Define: $f^a$

$f^a(⟨c_{10},\dots ,c_{1m}⟩,⟨c_{20},\dots ,c_{2m}⟩)$ compare to $f^a([l_1,u_2],[l_2,u_2])$

$⟹⟨c_{10}+c_{20},\dots ,c_{1m}+c_{2m}⟩$

$(0+{ϵ}_1,1+{ϵ}_2):f^a(⟨0,1,0⟩,⟨1,0,1⟩)=⟨1,1,1⟩$

output zonotope is the set $\{1+{ϵ}_1+{ϵ}_2|{ϵ}_1,{ϵ}_2\in [-1,1]\}$ which is the interval $[-1,3]$.

Affine Functions

$$f(x_1,\dots ,x_n)=\sum _j{a}_j{x}_j$$

where $a_j\in \mathbb{R}$.

$$f^a(⟨c_{1i}⟩,\dots ,⟨c_{ni}⟩)=⟨\sum _j{a}_j{c}_{j0},\dots ,\sum _j{a}_j{c}_{jm}⟩$$

$f(x,y)=3x+2y$
$f^a(⟨1,2,3⟩,⟨0,1,1⟩)=⟨f(1,0),f(2,1),f(3,1)⟩=⟨3,8,11⟩$
$3+8{ϵ}_1+11{ϵ}_2{ϵ}_1,{ϵ}_2\in [-1,1]$

Abstract Transformer of Activation Functions

For intervals domains we had: $rel{u}^a([l,u])=[relu(l),relu(u)]$
This formulation does not know how the inputs are related to their corresponding outputs. Geometrically the interval domain of ReLU approximates the function with a box as shown below: The breadth of the box depends on the lower bound of the interval, i.e., whether $l$ is positive or negative. Using zonotopes, the ReLU abstract transformer is build using a 1-dimensional zonotope $⟨c_i{⟩}_i$ as input

$$rel{u}^a(⟨c_i{⟩}_i)=\{\begin{array}{l}⟨c_i{⟩}_i\text{ for }l\ge 0\\ \\ ⟨0{⟩}_i\text{ for }u\le 0\\ \\ ?\text{ otherwise}\end{array}$$

If the lower bound is greater than zero, we return the input, if the upper bound is lesser than or equal to zero we return zero. Since zonotopes allows for relating input and output, we can shear the rectangles to form better approximations of ReLU.

The bottom face of the zonotope is $y=\lambda x$, for some slope $\lambda$. The top face is $y=\lambda x+u(1-\lambda )$. For $\lambda =0$, we get the base rectangle that is the interval domain for the ReLU, the steepness of the shear depends on the $\lambda$ parameter, which cannot be more than $u/(u-l)$ to ensure that the parallelogram covers the ReLU along the input range $[l,u]$. The distance between the top and bottom faces of the parallelogram is $u(1-\lambda )$, Thus the centre of the zonotope is at the point

$$\eta =\frac{u(1-\lambda )}{2}$$

With this information, we can now complete the definition of $rel{u}^a$ as follows:

$$rel{u}^a(⟨c_i{⟩}_i)=\{\begin{array}{l}⟨c_i{⟩}_i\text{ for }l\ge 0\\ \\ ⟨0{⟩}_i\text{ for }u\le 0\\ \\ ⟨\lambda c_1,...,\lambda c_m,0〉+⟨\eta ,0,0,...,\eta 〉\text{ otherwise}\end{array}$$

we have added a new generator, ${ϵ}_{m+1}$, in order to stretch the parallelogram in the vertical axis; its coefficient is $\eta$, which is half the height of the parallelogram. We also add the input zonotope scaled by $\lambda$ with coefficient $0$ for the new generator to ensure that we capture the relationship between the input and output.

Neural Polyhedron Abstraction

So far we have been approximating functions using a hyperrectangle with interval domains but it was non-relational, the zonotope domain allows us to approximate functions using a zonotope, e.g., a parallelogram and capture relations between different dimensions. The polyhedron domain is a more expressive abstract domain as it enables us to approximate functions using any arbitrary complex polyhedra. A polyhedron in ${\mathbb{R}}^n$ is a region made of straight faces. Convex polyhedra are shapes that have any two points of the shape completely contained in the shape and can be specified as a set of linear inequalities. This approximation is more precise for ReLU functions than those afforded by the interval and zonotope domains as the shapes used to approximated are not limited to hyperrectangle or parallelograms.

$$\begin{array}{rlrlrlrl}& \text{interval domain}& \to & & \text{zonotope abstract domain}& \to & & \text{polyhedra abstract domain}\\ & \text{abstracts a function}& & & \text{approximates a function using}& & & \text{approximates a function using}\\ & \text{using a hyperrectangle}& & & \text{zonotopes/paralleloprams}& & & \text{arbitrary convex polyhedra}\end{array}$$

Polyhedra are defined in a manner analogous to a zonotope abstractions, using a set of $m$ generator variables, ${ϵ}_1,\dots ,{ϵ}_m$ which are then bounded by a set of linear inequalities instead of being limited to the interval $[-1,1]$ as is the case for zonotopes.

A zonotope in ${\mathbb{R}}^n$

$$\left\{\left(c_{10}+\sum _{i=1}^m{c}_{1i}{ϵ}_i,\dots ,c_{n0}+\sum _{i=1}^m{c}_{ni}{ϵ}_i\right)|F({ϵ}_1,\dots ,{ϵ}_m)\right\}$$

where $F({ϵ}_1,\dots ,{ϵ}_m)$ is a Boolean statement that evaluates to $true$ iff all of its arguments are in $[-1,1]$. For a polyhedron $F$ is defined as a set (conjunction) of linear inequalities over the generator variables, e.g.,

$$F({ϵ}_1,{ϵ}_2)\equiv 0\le {ϵ}_1\le 5\wedge {ϵ}_1={ϵ}_2$$

$F$ defines a bounded polyhedron over the generator variables, giving a lower and upper bound for each generator, e.g., ${ϵ}_1\le 0$ is not allowed, because it does not enforce a lower bound on ${ϵ}_1$. In the 1-dimension, a polyhedron is simply an bounded interval.

To find the upper and lower bounds of a polyhedron, we need to solve a linear program which takes polynomial time in the number of variables and constraints.

To compute the lower bound of the $j$-th dimension, we solve for the following:

$$\text{min }{c}_{j0}+\sum _{i=1}^m{c}_{ij}{ϵ}_i$$

Similarly to calculate the upper bound of $j$-th dimension, we just take the $\text{max}$ instead of minimizing the generator terms.

A given polyhedra in ${\mathbb{R}}^n$

$$\left\{\left(c_{10}+\sum _{i=1}^m{c}_{1i}{ϵ}_i,\dots ,c_{n0}+\sum _{i=1}^m{c}_{ni}{ϵ}_i\right)|F({ϵ}_1,\dots ,{ϵ}_m)\right\}$$

will be abbreviated as:

$$(⟨c_{1i}{⟩}_i,\dots ,⟨c_{ni}{⟩}_i,F)$$

Example: Given polyhedron: $\{({ϵ}_1,{ϵ}_2)|F({ϵ}_1,{ϵ}_2)\}$

$F\equiv (0\le {ϵ}_1\le 1)\wedge ({ϵ}_2\le {ϵ}_1)\wedge ({ϵ}_2\ge 0)$

From the Boolean function we can understand that a bounded x-y space is being defined. For ease of visualisation, we can replace ${ϵ}_1$ and ${ϵ}_2$ with $x$ and $y$ which gives us a set pf equations that $x$ and $y$ must satisfy:

$$\begin{array}{rl}& 0\le x\le 1\\ & y\le x\\ & y\ge 0\end{array}$$

This set of equations clearly defines a LP problem in the shape of a triangle situated in the first quadrant of the x-y plane. Thus the polyhedron comprises of all of the points inside this triangle.

$(⟨0,1,0⟩,⟨0,0,1⟩,F)$

We need to solve for the following four equations to get the upper and lower bound on the generators

$$\begin{array}{rlrlrlrlr}& \text{max }{ϵ}_1& & \text{min }{ϵ}_1& & \text{max }{ϵ}_2& & \text{min }{ϵ}_2& \\ & \text{s.t. }F& & \text{s.t. }F& & \text{s.t. }F& & \text{s.t. }F& \end{array}$$

Thus the valid interval domains of the generators are: ${ϵ}_1\in [0,1]$ and ${ϵ}_2\in [0,1]$.

Abstract Transformers for Polyhedra

Affine Functions

For an affine function $f(x_1,\dots ,x_n)={\sum }_j{a}_j{x}_j{a}_j\in \mathbb{R}$ We will have the following abstract transformer:

$$f^a(⟨c_{1i}⟩,⟨c_{ni}⟩,F)=\left(⟨\sum _j{a}_j{c}_{j0},\dots ,\sum _j{a}_j{c}_{jm}⟩,F\right)$$

the set of linear equalities does not change between input and output of the function.

Example:

$f(x,y)=3x+2y$

$f^a(⟨1,2,3⟩,⟨0,1,1⟩,F)=(⟨3,8,11⟩,F)$

The abstract transformer for ReLU generated in polyhedron domain:

$rel{u}^a$ in convex polyhedra

Top face: $y=\frac{u(x-l)}{u-l}$ from the equation of a straight line, where $x_1=l$, $x_2=u$ and $y_1=relu(l)$, $y_2=relu(u)$.

We need to compute the shaded area which is bounded by $y=0$ from below, $y=x$ from the right, and $y=\frac{u(x-l)}{u-l}$ from above. We define $rel{u}^a$ as

$$rel{u}^a(⟨c_i{⟩}_i,F)=(⟨\underset{m}{\underbrace{0,0,\dots ,0}},1⟩,F^{\prime })$$

where

$$F^{\prime }\equiv F\wedge \left({ϵ}_{m+1}\le \frac{u(⟨c_i⟩-l)}{(u-l)}\right)\wedge ({ϵ}_{m+1}\ge 0)\wedge ({ϵ}_{m+1}\ge ⟨c_i⟩)$$

• $l$ and $u$ are the lower and upper bounds of the input polyhedron that can be computed using linear programming.
• $〈c_i{〉}_i$ is being used to denote the full term $c_0+{\sum }_{i=1}^m{c}_i{ϵ}_i$.
• And a new generator, ${ϵ}_{m+1}$ has been added, the new set of constraints $F^{\prime }$ relates this new generator to the input, effectively defining the shaded region.

Condensing the constraints into two dimensions:

$$rel{u}^a(⟨0,1⟩,-1\le {ϵ}_1\le 1)=(⟨0,0,1⟩,F^{\prime })$$

where ${ϵ}_1$ determines $x$ and ${ϵ}_2$ determines $y$ thus,

$$F^{\prime }\equiv \left(-1\le {ϵ}_1\le 1\right)\wedge \left({ϵ}_2\le \frac{{ϵ}_1+1}{2}\right)\wedge ({ϵ}_2\ge 0)\wedge ({ϵ}_2\le {ϵ}_1)$$

Abstractly Interpreting Neural Networks

$G=(V,E)f_G:{\mathbb{R}}^n\to {\mathbb{R}}^m$

$$f_G^a:\text{n-dimensional polyhedron}\to \text{m-dimensional polyhedron}$$

$f_G^a(⟨c_{1j}⟩,\dots ,⟨c_{nj}⟩,F)$

• For every input node $v_i$, we have $ou{t}^a(v_i)=(⟨c_{ij}{⟩}_j,F)$

• For every non-input node $v$ we have $ou{t}^a(v)=f_v^a(p_1,\dots ,p_k,{\bigwedge }_{i=1}^n{F}_k)$ where $f_v^a$ is the abstract transformer of $f_v$, $v$ has incoming edges $(v_1,v),\dots ,(v_k,v)$ and $ou{t}^a(v_i)=(p_i,F_i)$

• output of $f_G^a$ is the $m$-dimensional polyhedron $\left(p_1,\dots ,p_m,{\bigwedge }_{i=1}^m{F}_i\right)$, where $v_1,\dots ,v_m$ are the output nodes and $ou{t}^a(v_i)=(p_i,F_i)$

Some abstract transformers for activation functions add new generators, we assume all of them were already in the polyhedron but with coefficients set to $0$, they get non-zero coefficients only in the output of activation functions.

Abstract Interpretation based Verification

For verification using the Hoare triplet to work, we require the a sound representation of the set of val of $x$ in the abstract domain. Along with the abstract representation of the neural network $f$ on all values of $x$ that results in an over-approximation of the values of $x$. We also need to check that all values of $r$ satisfy the postcondition.

The generic precondition can be defined as:

$$\begin{array}{r}\{\Vert x-c{\Vert }_p\le ϵ\}\\ r\leftarrow f(x)\\ \{\text{class}(r)=y\}\end{array}$$

Here for the sake of generality, we do not specify a particular norm, instead we define an $l_p$-norm. The main norms that we will be considering in this case are the $l_2$-norm and the $l_{\infty }$-norm both of which are distance metrics but serve very different purposes.

$$l_p\text{ norm: }\Vert z{\Vert }_p\{\begin{array}{l}{l}_2\text{ norm: }\Vert z{\Vert }_2=\sqrt{{\sum }_i|z_i{|}^2}\\ \\ l_{\infty }\text{ norm: }\Vert z{\Vert }_{\infty }=\underset{i}{\max }|z_i|\end{array}$$

$l_2$-norm is the length of the straight line between two images in ${\mathbb{R}}^n$ and $l_{\infty }$-norm is the largest discrepancy between two corresponding pixels. They have their own usability based on what kinds of deviations we wish to categorize. If there is a great deal of variability restricted to some locality of the image we can use $l_2$-norm as it allows a small number of pixels to significantly differ in brightness. If the noise is random and spread out through the entire image, $l_{\infty }$-norm is more suitable as it bounds the maximum discrepancy a corresponding pixels in the two images can have.

Abstracting the Precondition

We have the precondition as: $\{x|\Vert x-c{\Vert }_{\infty }\le ϵ\}$

Interval: $I=([c_1-ϵ,c_1+ϵ],\dots ,[c_n-ϵ,c_n+ϵ])$ $\Vert \cdot {\Vert }_{\infty }$ allows us to take elements of $c$ and change it by $ϵ$ independently of other dimensions. Define abstract transformer $f^a(I)=([l_1,u_1],[l_2,u_2],\dots ,[l_m,u_m])$

$I^{\prime }=([l_1,u_1],[l_2,u_2],\dots ,[l_m,u_m])$ represents all possible values of $r$ and more.

We have to prove that $\forall r\in I^{\prime }$, $class(r)=y$ We see that if $l_y>u_i\forall i\ne y$ then, $\forall r\in I^{\prime },class(r)=y$

If $y$-th interval is larger than all others, then we know that the classification is always $y$.

Note: if $l_y\le u_i$ for some $i\ne y$, then we cannot disprove this property, so this is a one-sided check. In simpler words, since we will be receiving a range of probabilities corresponding to each of the class indices. We will be choosing the class whose lower bound probability ($l_y$) is higher than all other upper bounds of all other class ($u_i$). We will be picking the non-overlapping bound with highest lower bound. If there is any overlap between the range of probabilities, we cannot be entirely sure as to which class the resulting probability belongs to inside of the overlapping region. Thus making this is a one-sided check. Example:

$f^a(I)=I^{\prime }=([0.1,0.2],[0.3,0.4])$ $class(r)=2\forall r\in I^{\prime }$

where $I^{\prime }=([0.1,0.2],[0.15,0.4])$

two intervals overlap in the 0.15 to 0.2 region, this means that we cannot conclusively say that $class(r)=2\forall r\in I^{\prime }$, so verification fails. $\exists r\in I^{\prime }$ that can belong to both $[0.1,0.2]$ and $[0.15,0.4]$.

Verifying Robustness with Zonotopes

Checking $l_{\infty }$ robustness property using zonotopes. Since precondition is hyperrectangular it can be precisely represented:

$$f^a(z)=z^{\prime }$$

we want to ensure that the dimension $y$ is greater than all others, problem is akin to checking if a 1-D zonotope is always > 0.

$$z^{\prime }=(⟨c_{1i}⟩,\dots ,⟨c_{mi}⟩)$$

To check that dimension $y$ is greater than dimension $j$, we check if the lower bound of the 1D zonotope $⟨c_{yi}⟩-⟨c_{ji}⟩$ is > 0 or not

Example: $z^{\prime }=(2+{ϵ}_1,4+{ϵ}_1+{ϵ}_2)$ For this region, we can graphically see that $y>x$ for any point $(x,y)$ To check $y>x$ mechanically, we subtract the $x$-dimension from the $y$-dimension.

$$(4+{ϵ}_1+{ϵ}_2)-(2+{ϵ}_1)=2+{ϵ}_1$$

The resulting 1D zonotope $(2+{ϵ}_1)$ denotes the interval $[1,3]$ which is greater than zero.

So, for multiple $m$-generators and $n$-dimensions

$$(⟨c_0+c_1{ϵ}_1+\cdots +c_m{ϵ}_m{⟩}_1,⟨c_0+c_1{ϵ}_1+\cdots +c_m{ϵ}_m{⟩}_2,\dots ,⟨c_0+c_1{ϵ}_1+\cdots +c_m{ϵ}_m{⟩}_n)$$

for some $y$-dimension check $⟨c_y{⟩}_y-⟨c_i{⟩}_i$ for every $i\ne y$ and check whether $>0$

(we will be getting a set of ranges, if every range > 0 then $y$ is the class)

Verifying Robustness with Polyhedra

Analogous to zonotopes but requires invoking a linear program solver.

We represent the interval as a hyperrectangle polyhedron $Y$. Then we evaluate $f^a(Y)$ resulting in a polyhedron

$$Y^{\prime }=(⟨c_{1i}⟩,\dots ,⟨c_{mi}⟩,F)$$

To check if dimension $y$ is greater than dimension $j$, we ask a linear program-solver if the following constraints are satisfiable

$$F\wedge ⟨c_{yi}⟩>⟨c_{ji}⟩$$

Robustness in $l_2$-norm

$\{x|\Vert x-c{\Vert }_2\le ϵ\}$ : this defines a unit circle around $(0,0)$.

let $c=(0),ϵ=1$

Cannot be precisely represented in the interval domain, best we can do is $([-1,1],[-1,1])$ with polyhedra, we can define polyhedra with more and more faces to more accurately approximate a circle, but there is a precision-scalability trade-off that comes into the picture.

Robustness in NLP

Synonyms of words should not confuse Neural Networks. Complete set of synonyms: $S_w$ ; each word is $w$ $i$-th element of a vector is the $i$-th token/word embedding of the sentence

Correctness Property:

$$\begin{array}{r}\{x_i\in S_{c_i}\forall i\}\\ r\leftarrow f(x)\\ \{class(r)=y\}\end{array}$$

all vectors $x$ that are like $c$ but where some words are replaced by synonyms. Set of possible vectors $x$ is finite but exponential in length of input sentences

$$([\text{min }{S}_{c_1},\text{max }{S}_{c_1}],\dots ,[\text{min }{S}_{c_n},\text{max }{S}_{c_n}])$$

Abstract Training of Neural Networks

Training Data: $\{(x_1,y_1),\dots ,(x_m,y_m)\}$ $x_i\in {\mathbb{R}}^n$ binary label (classification model): $y_i\in \{0,1\}$

We assume we have a family of functions represented as a parameterized function: $f_{\theta }$ Where $\theta$ is the vector of weights. Search the space of $\theta$ and find the best values of $\theta$.

Family of affine functions: $f_{\theta }(x)={\theta }_1+{\theta }_2{x}_1+{\theta }_3{x}_2$

Solve optimization problem: $\text{arg }\underset{\theta }{\min }\frac{1}{m}{\sum }_{i=1}^{m}1\left[f_{\theta }(x_i)=y_i\right]$

$$1[b]\{\begin{array}{l}1\text{ if b is True}\\ \\ 0\text{ if b is False}\end{array}$$

This function is difficult to resolve as Boolean is non-differential so we use MSE:

$$\text{arg }\underset{\theta }{\min }\frac{1}{m}\sum _{i=1}^m{\left(f_{\theta }(x_i)-y_i\right)}^2$$

family of functions $f_{\theta }$ represented as neural net graph $G_{\theta }$, where every node $v$‘s function $f_v$ may be parameterized by $\theta$. Formally we solve:

$$\text{arg }\underset{\theta }{\min }\frac{1}{m}\sum _{i=1}^{m}L(\theta ,x_i,y_i)$$

Loss function $L$ can also be represented as a neural network. $L$ is represented as an extension of the graph $G_{\theta }$ by adding a node at the very end that computes the loss function.

$f_{\theta }:{\mathbb{R}}^n\to \mathbb{R}$

flowchart LR

v1[v₁] & v2[v₂] --> A((...)) & B((...))
A & B--> vo[vₒ]

subgraph "intermediate nodes"
A
B
end

We can construct a graph for the loss function $L(\theta ,x,y)$ by adding am input node $v_y$ for the label $y$ and creating a new output node $v_L$ that compares the output of $f_{\theta }$ (the node $v_{\theta }$) with $y$.

flowchart LR

v1[v₁] & v2[v₂] -.- vo((vₒ)) --> vL[v_L]
vy[v_y] --> vL

input node $v_y$ takes in the label $y$ and $f_{v_L}$ encodes the loss function, e.g., MSE

$$\left(\frac{\partial g}{\partial {\theta }_1},\dots ,\frac{\partial g}{\partial {\theta }_n}\right)\overrightarrow{\text{ for a particular }{\theta }^0}\left(\frac{\partial g({\theta }^0)}{\partial {\theta }_1},\dots ,\frac{\partial g({\theta }^0)}{\partial {\theta }_n}\right)$$

Gradient Descent:

• Start with $j=0$ and a random value $\theta$, called ${\theta }^0$
• Set ${\theta }^{j+1}$ to ${\theta }^j-\eta ((\nabla g)({\theta }^i))$
• Set $j$ to $j+1$ and repeat $\eta >0$ learning rate

$$\frac{1}{m}\sum _{i=1}^m\nabla L(\theta ,x_i,y_i)$$

Set ${\theta }^{j+1}$ to ${\theta }^j-\frac{\eta }{m}{\sum }_{i=1}^m\nabla L({\theta }^j,x_i,y_i)$

SGD/mini-batch gradient descent:

• Start with $j=0$ and a random value $\theta$, called ${\theta }^0$
• Divide the dataset into a random set of $k$ batches: $B_1,B_2,\dots ,B_k$
• For $i$ from 1 to $k$:
  • Set ${\theta }^{j+1}$ to ${\theta }^j-\frac{\eta }{m}{\sum }_{(x,y)\in B_i}^m\nabla L({\theta }^j,x_i,y_i)$
  • Set $j$ to $j+1$
• reiterate

Size of batches $k$ is dependent on how much data can be put into the GPU at any one point.

These optimization algorithms do not provide any robustness to the networks, as we are only minimizing loss over average prediction. Not robust to perturbations, even if they are, abstract interpretation fails to provide a proof due to its overapproximate nature. So we would like to train neural networks in such a way that they are friendly to abstract interpretation.

Redefining Optimization Objective for Robustness

(Robustness Optimization Objective)

For every $(x,y)$ in our dataset, we want the neural net to predict $y$ on all images $z$ such that $\Vert x-z{\Vert }_{\infty }\le ϵ$. This set can be characterized as

$$R(x)=\{z|\Vert x-z{\Vert }_{\infty }\le ϵ\}$$

New Optimization Objective

$$\text{arg }\underset{\theta }{\min }\frac{1}{m}\sum _{i=1}^m\underset{z\in R(x_i)}{\max }L(\theta ,x_i,y_i)$$

Instead of minimizing for the loss of $(x_i,y_i)$, we minimize loss for the worst-case perturbation of $x_i$ from the set $R(x_i)$. This is known as robust optimization problem. Training the neural net using such an objective is known as adversarial training (very similar in appearance to minimax or other game playing techniques).

Solving Robust Optimization via Abstract Interpretation

$R(x)$ can be defined in interval domain precisely as it represents a set of images within an $l_{\infty }$ bound. We can overapproximate the inner maximization by abstractly interpreting $L$ on the entire set $R(x_i)$. By virtue of the soundness of the abstract transformer $L^a$, we know:

$$\left(\underset{z\in R(x_i)}{\max }L(\theta ,x_i,y_i)\right)\le u$$

where $L^a(\theta ,R(x_i),y_i)=[l,u]$ Therefore we can overapproximate the inner maximization by abstractly interpreting the loss function on the set $R(x_i)$ and taking the upper bound.

Thus the robust optimization objective can be rewritten as:

$$\text{arg }\underset{\theta }{\min }\frac{1}{m}\sum _{i=1}^m\text{upper bound of }{L}^a(\theta ,R(x_i),y_i)$$

instead of treating $L^a$ as an abstract transformer, we can treat it as a function taking a vector of inputs and returning upper and lower bounds, this is called flattening the abstract transformer.

Example: $relu(x)=max(0,x)$

$rel{u}^a([l,u])=[max(0,l),max(0,u)]$

flattening the $rel{u}^a$ to $rel{u}^{af}:{\mathbb{R}}^2\to {\mathbb{R}}^2$

$rel{u}^a([l,u])=(max(0,l),max(0,u))$ : returns a tuple of values instead of an interval.

$$\text{arg }\underset{\theta }{\min }\frac{1}{m}\sum _{i=1}^m{L}_u^{af}(\theta ,l_{i1},u_{i1},\dots ,l_{in},u_{in},y_i)$$

Where $L_u^{af}$ is the only upper bound of $L^{af}$ output

$$R(x_i)=([l_{i1},u_{i1}],\dots ,[l_{in},u_{in}])$$

SGD can optimize for such objectives because all of the abstract transformers of the interval domain that are of interest for neural nets are differentiable almost everywhere. Some can be adapted into zonotopes also.

Example: $f:\mathbb{R}\to \mathbb{R}$

$$\begin{array}{rlrlrlr}& f^a:& & \text{1D zonotope with m-generators}& \to & & \text{1D zonotope with m-generators}\\ & & & ⟨c_0,c_1,c_2,c_3,c_4,c_5,\dots ,c_m⟩& & & ⟨c_0^{\prime },c_1^{\prime },c_2^{\prime },c_3^{\prime },c_4^{\prime },c_5^{\prime },\dots ,c_m^{\prime }⟩\end{array}$$

$f^{af}:{\mathbb{R}}^{m+1}\to {\mathbb{R}}^{m+1}$

Flattening does not work for polyhedra domain, because it invokes a blackbox linear programming solver for activation functions which is not differentiable.

Neural Networks trained with abstract interpretation are:

• more robust to perturbation attacks
• verifiably robust using abstract interpretation A Neural Network could satisfy a correctness property but we may not be able to verify the neural net using abstract domain. By incorporating abstract interpretation right into the training we guide SGD towards neural nets that are more amenable to verification.

Note: $l_p$ robustness properties are closely related to the notion of Lipschitz continuity. For instance, a network $f:{\mathbb{R}}^n\to {\mathbb{R}}^m$ is K-Lipschitz under the $l_2$-norm if

$$\Vert f(x)-f(y){\Vert }_2\le K\Vert x-y{\Vert }_2$$

The smallest K satisfying the above is called the Lipschitz constant of $f$. If we can bound K, then we can prove $l_2$-robustness of $f$.