Mixed Integer Linear Programs for Adversarial Robustness

\[\newcommand{\kl}[1]{\mathopen{}\left( #1 \right)\mathclose{}} \newcommand{\ekl}[1]{\mathopen{}\left[ #1 \right]\mathclose{}} \newcommand{\skl}[1]{\mathopen{}\left\{ #1 \right\}\mathclose{}} \newcommand{\bkl}[1]{\mathopen{}\left| #1 \right|\mathclose{}} \newcommand{\nkl}[1]{\mathopen{}\left\| #1 \right\|\mathclose{}} \newcommand{\bfa}{\mathbf{a}} \newcommand{\bfb}{\mathbf{b}} \newcommand{\bfc}{\mathbf{c}} \newcommand{\bfd}{\mathbf{d}} \newcommand{\bfe}{\mathbf{e}} \newcommand{\bff}{\mathbf{f}} \newcommand{\bfg}{\mathbf{g}} \newcommand{\bfh}{\mathbf{h}} \newcommand{\bfi}{\mathbf{i}} \newcommand{\bfj}{\mathbf{j}} \newcommand{\bfk}{\mathbf{k}} \newcommand{\bfl}{\mathbf{l}} \newcommand{\bfm}{\mathbf{m}} \newcommand{\bfn}{\mathbf{n}} \newcommand{\bfo}{\mathbf{o}} \newcommand{\bfp}{\mathbf{p}} \newcommand{\bfq}{\mathbf{q}} \newcommand{\bfr}{\mathbf{r}} \newcommand{\bfs}{\mathbf{s}} \newcommand{\bft}{\mathbf{t}} \newcommand{\bfu}{\mathbf{u}} \newcommand{\bfv}{\mathbf{v}} \newcommand{\bfw}{\mathbf{w}} \newcommand{\bfx}{\mathbf{x}} \newcommand{\bfy}{\mathbf{y}} \newcommand{\bfz}{\mathbf{z}} \newcommand{\bfA}{\mathbf{A}} \newcommand{\bfB}{\mathbf{B}} \newcommand{\bfC}{\mathbf{C}} \newcommand{\bfD}{\mathbf{D}} \newcommand{\bfE}{\mathbf{E}} \newcommand{\bfF}{\mathbf{F}} \newcommand{\bfG}{\mathbf{G}} \newcommand{\bfH}{\mathbf{H}} \newcommand{\bfI}{\mathbf{I}} \newcommand{\bfJ}{\mathbf{J}} \newcommand{\bfK}{\mathbf{K}} \newcommand{\bfL}{\mathbf{L}} \newcommand{\bfM}{\mathbf{M}} \newcommand{\bfN}{\mathbf{N}} \newcommand{\bfO}{\mathbf{O}} \newcommand{\bfP}{\mathbf{P}} \newcommand{\bfQ}{\mathbf{Q}} \newcommand{\bfR}{\mathbf{R}} \newcommand{\bfS}{\mathbf{S}} \newcommand{\bfT}{\mathbf{T}} \newcommand{\bfU}{\mathbf{U}} \newcommand{\bfV}{\mathbf{V}} \newcommand{\bfW}{\mathbf{W}} \newcommand{\bfX}{\mathbf{X}} \newcommand{\bfY}{\mathbf{Y}} \newcommand{\bfZ}{\mathbf{Z}} \newcommand{\bfone}{\mathbf{1}} \newcommand{\bfzero}{\mathbf{0}} \newcommand{\E}{\mathbb{E}} \newcommand{\R}{\mathbb{R}} \renewcommand{\P}{\mathbb{P}} \newcommand{\bfmu}{\bm{\mu}} \newcommand{\bfsigma}{\bm{\sigma}} \newcommand{\bfdelta}{\boldsymbol{\delta}} \newcommand{\bfSigma}{\bm{\Sigma}} \newcommand{\bfLambda}{\bm{\Lambda}} \newcommand{\bfeta}{\bm{\eta}} \newcommand{\bftheta}{\bm{\theta}} \newcommand{\CA}{\mathcal{A}} \newcommand{\CB}{\mathcal{B}} \newcommand{\CC}{\mathcal{C}} \newcommand{\CD}{\mathcal{D}} \newcommand{\CE}{\mathcal{E}} \newcommand{\CF}{\mathcal{F}} \newcommand{\CG}{\mathcal{G}} \newcommand{\CH}{\mathcal{H}} \newcommand{\CI}{\mathcal{I}} \newcommand{\CJ}{\mathcal{J}} \newcommand{\CK}{\mathcal{K}} \newcommand{\CL}{\mathcal{L}} \newcommand{\CM}{\mathcal{M}} \newcommand{\CN}{\mathcal{N}} \newcommand{\CO}{\mathcal{O}} \newcommand{\CP}{\mathcal{P}} \newcommand{\CQ}{\mathcal{Q}} \newcommand{\CR}{\mathcal{R}} \newcommand{\CS}{\mathcal{S}} \newcommand{\CT}{\mathcal{T}} \newcommand{\CU}{\mathcal{U}} \newcommand{\CV}{\mathcal{V}} \newcommand{\CW}{\mathcal{W}} \newcommand{\CX}{\mathcal{X}} \newcommand{\CY}{\mathcal{Y}} \newcommand{\CZ}{\mathcal{Z}} \newcommand{\frA}{\mathfrak{A}} \newcommand{\frB}{\mathfrak{B}} \newcommand{\frC}{\mathfrak{C}} \newcommand{\frD}{\mathfrak{D}} \newcommand{\frE}{\mathfrak{E}} \newcommand{\frF}{\mathfrak{F}} \newcommand{\frG}{\mathfrak{G}} \newcommand{\frH}{\mathfrak{H}} \newcommand{\frI}{\mathfrak{I}} \newcommand{\frJ}{\mathfrak{J}} \newcommand{\frK}{\mathfrak{K}} \newcommand{\frL}{\mathfrak{L}} \newcommand{\frM}{\mathfrak{M}} \newcommand{\frN}{\mathfrak{N}} \newcommand{\frO}{\mathfrak{O}} \newcommand{\frP}{\mathfrak{P}} \newcommand{\frQ}{\mathfrak{Q}} \newcommand{\frR}{\mathfrak{R}} \newcommand{\frS}{\mathfrak{S}} \newcommand{\frT}{\mathfrak{T}} \newcommand{\frU}{\mathfrak{U}} \newcommand{\frV}{\mathfrak{V}} \newcommand{\frW}{\mathfrak{W}} \newcommand{\frX}{\mathfrak{X}} \newcommand{\frY}{\mathfrak{Y}} \newcommand{\frZ}{\mathfrak{Z}} \newcommand{\CNP}{\mathcal{NP}} \newcommand{\CPP}{\mathcal{PP}} \newcommand{\SP}{\mathsf{P}} \newcommand{\SPP}{\mathsf{PP}} \newcommand{\SSP}{\mathsf{\#P}} \newcommand{\SNP}{\mathsf{NP}} \newcommand{\SBPP}{\mathsf{BPP}} \newcommand{\ScoNP}{\mathsf{coNP}} \newcommand{\bbone}{\mathbbm{1}} \newcommand{\ord}{\mathrm{ord}} \newcommand{\odr}{\vee} \newcommand{\und}{\wedge} \newcommand{\Odr}{\bigvee} \newcommand{\Und}{\bigwedge} \newcommand{\xor}{\oplus} \newcommand{\Xor}{\bigoplus} \newcommand{\bmat}[1]{\begin{bmatrix} #1 \end{bmatrix}} \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator*{\argmax}{argmax}\]

The Problem:

The problem we try to solve is

\[\max_{\bfx \in C} f(\bfx),\]

where $f$ is given as a neural network and $C$ is some convex region. This problem arises from questions of adversarial robustness, where $ C = \skl{\bfx\in \R^n:~ \nkl{\bfx - \bfx_{\text{sample}} } \leq \epsilon}, $ formal verification of networks and interpretability.

Generally, this is a $\SNP$-hard problem and any general solution strategy will take exponentially long in the worst case. However, real-world trained neural networks do not necessarily represent the worst-case. Thus a heuristic algorithm like Gradient Descent works reasonably well in solving this problem.

However, gradient descent does not give us any bounds on how close we are to the optimal $\bfx$. We would like to design an algorithm that combines formal guarantees with speedy solutions for most real-world instances.

Using MIP Solvers

A ReLU neural network can be rewritten as a Mixed Integer Linear Program when there exists a general bound on the size of the variables in each layber, see this paper for instructions.

The general idea is to introduce binary variables that signify on the sign of the pre-Rely neuron and thus which linear region of the ReLU is applied. When all these sign variables are fixed the neural network becomes a linear problem.

Of course, searching over all possible assignments takes exponentially long in the number of ReLU-activated neurons. Solving the whole network at once might therefore be only feasible for small problems.

The Constraint Propagation Algorithm (CPA)

This is where the new constrain propagation algorithm comes in.

Preparation:

Get a rough first bound on all the post-Relu activations $x^{(l)}_j$ in every layer $l$ by using the inequality

\[0 \leq x^{(l)}_j \leq c^{(l)}_j\]

with

\[c^{(l)}_j = \bfw^{(l),+}_j \bfc^{(l)} + b^{(l)}_j\]

and

\[c^{(1)}_j = \sum_i \bfw^{(1),+}_{ij} \max(x_i) + \bfw^{(1),-}_{ij} \min(x_i),\]

where $\bfw^{(l),+}$ and $\bfw^{(l),-}$ are the positive and negative part of the weight matrix in layer $l$ and $\max(x_i)$ and $\min(x_i)$ are the upper and lower bound on the input variables.

This will give us a reasonable a priori bound on all the activations. This bound is then used to transform the network into a MIP.

The idea of the CPA is to iteratively make these bounds tighter in the relevant direction.

The algorithm:

The algorithm loops over alternating forward and a backward pass.

Forward Pass:

Propagate a feasible candidate vector from input to output layer and note the resulting candidate for every layer. In the beginning, this vector can be randomly chosen. We can update a lower bound on the target value in the output layer with each forward pass.

Backward Pass:

Now, we start at the output layer and go back to the input.

For layer l we assume to have:
- A target vector $\bft$
- A candidate vector $\bfk$
- A convex bounding region $C$ for the neurons in layer $l-1$
First, try to solve

\[\bfx \in C \quad \text{s.t.}\quad \bft = \text{ReLU}(W^{(l)T} \bfx + \bfb^{(l)})\]

a) Feasible case: Proceed to layer $l-1$ with new target vector $\bfx$.

b) Infeasible case:

Determine active direction $\bfd = \bft - \bfk$.
Solve linear problem

\[\max_{\bfx \in C} \bfd^T \text{ReLU}(W^{(l)T} \bfx + \bfb^{(l)})\]

which gives solution $\bfx^\ast$, value $v^\ast$ and a set of active constraints.

We can thus add

\[\bfd^{T} \bfx \leq v^*\]

as a new constraint for layer $l$.

Proceed to layer $l-1$ with new target $\bfx^*$

We start at the output layer $L$ with the (1-dimensional) target $c^L$, and repeat the procedure once we reached the input layer. Note that the constraints on every layer get tighter over time, which leads to increasingly realistic target vectors.

The forward pass increases the lower bound on the maximal output value, while the backward pass decreases the upper bound. We are done when both converge or are sufficiently close.

Figure 1. Illustation for a layer of width 2. The known constraints result in the light green convex area. The dark green arrows mark the linear problem for each ReLU region. The point $a^\ast$ is the solution to the resulting maximisation problem. If this target cannot be reached in the next lower layer, the active direction $\bfd = \bfa^* - \bfk$ is optimised instead, which leads to the next maximisation problem.}

How to determine the active direction?

This is where a bit of experimentation comes in what works best. My first best guess is to chose as the difference between the candidate and the target. Another possibility is to consider the active constraints for the target vector and choose a vector in a way that maximally reduces the size of the convex region in that corner.

How to get started?

I think a good start is a neural network of width 2, where the constraints on every level can be easily visualised and more intuition can be gained. An easy beginner task might be a neural network of a single input and with depth 5 that fits some analytical function, like $\sin(x)$ over some compact domain, like $[0,2\pi]$.

Combination with ResNet Architecture:

The speed of each single MIP-solving step depends on the number of ReLU neurons. Neurons with a simple linear activation do not make the MIP problem harder. The ResNet architecure mixes neurons with ReLU acitvations and linear activations. We might reach a sweet spot for training performance and solver speed by trading the depth of the network against the number of ReLU units per layer.

Research Questions:

Is the layer-wise strategy faster than just solving the whole DNN-MIP at once?
What are the most complex datasets that can be solved in reasonable time?
How should we choose the active direction?
What is the optimal ResNet trade-off?
Prove that the upper and lower bounds are always correct (I think this is straightforward).