Safety Verification and Robustness Analysis of Neural Networks via Quadratic Constraints and Semidefinite Programming

TL;DR

This paper introduces a semidefinite programming framework using quadratic constraints to certify neural network safety and robustness against input uncertainties and adversarial attacks, applicable to general activation functions.

Contribution

It presents a novel SDP-based method that abstracts activation functions with quadratic constraints for safety verification of neural networks.

Findings

Effective in bounding neural network outputs under input uncertainties.

Balances conservatism and computational efficiency.

Applicable to various neural network architectures and activation functions.

Abstract

Certifying the safety or robustness of neural networks against input uncertainties and adversarial attacks is an emerging challenge in the area of safe machine learning and control. To provide such a guarantee, one must be able to bound the output of neural networks when their input changes within a bounded set. In this paper, we propose a semidefinite programming (SDP) framework to address this problem for feed-forward neural networks with general activation functions and input uncertainty sets. Our main idea is to abstract various properties of activation functions (e.g., monotonicity, bounded slope, bounded values, and repetition across layers) with the formalism of quadratic constraints. We then analyze the safety properties of the abstracted network via the S-procedure and semidefinite programming. Our framework spans the trade-off between conservatism and computational efficiency and applies to problems beyond safety verification. We evaluate the performance of our approach via numerical problem instances of various sizes.

Tables

Table 1. TABLE I : Average values (over 100 runs) of different upper bounds for the problem sup x ∈ 𝒳 f ​ ( x ) subscript supremum 𝑥 𝒳 𝑓 𝑥 \sup_{x\in\mathcal{X}}f(x) with 𝒳 = ‖ x − x ⋆ ‖ ∞ ≤ ϵ 𝒳 subscript norm 𝑥 subscript 𝑥 ⋆ italic-ϵ \mathcal{X}=\|x-x_{\star}\|_{\infty}\leq\epsilon , x ⋆ = 1 n x subscript 𝑥 ⋆ subscript 1 subscript 𝑛 𝑥 x_{\star}=\mathrm{1}_{n_{x}} and ϵ = 0.2 italic-ϵ 0.2 \epsilon=0.2 . The neural network f 𝑓 f has n x = 10 subscript 𝑛 𝑥 10 n_{x}=10 inputs, n f = 1 subscript 𝑛 𝑓 1 n_{f}=1 output and ℓ ∈ { 1 , ⋯ , 5 } ℓ 1 ⋯ 5 \ell\in\{1,\cdots,5\} hidden layers.

	Bounds				Running Time (Sec)
$\ell$	$\mathrm{MILP}$	$\mathrm{DeepSDP}$	$\mathrm{SDR}$	$\mathrm{LP}$	$\mathrm{MILP}$	$\mathrm{DeepSDP}$	$\mathrm{SDR}$	$\mathrm{LP}$
$1$	$1.07$	$1.12$	$1.13$	$1.81$	$0.04$	$0.82$	$0.55$
$2$	$2.04$	$2.52$	$2.74$	$7.62$	$25.96$	$8.26$	$4.71$
$3$	-	$11.08$	$12.21$	$50.60$	-	$34.18$	$31.20$
$4$	-	$47.74$	$54.15$	$368.65$	-	$78.95$	$94.74$
$5$	-	$218.8$	$266.3$	$3004.9$	-	$164.63$	$207.77$

Equations

$$\frac{\alpha \beta }{\alpha +\beta }\Vert y-x{\Vert }_2^2+\frac{1}{\alpha +\beta }\Vert \nabla g(y)-\nabla g(x){\Vert }_2^2$$

$$\le (\nabla g(y)-\nabla g(x){)}^{\top }(y-x),$$

$$\mathcal{Y}=f(\mathcal{X}):=\{y\in {\mathbb{R}}^{n_f}\mid y=f(x),x\in \mathcal{X}\}.$$

$$x^0$$

$$x^{k+1}$$

$$f(x)$$

$$\varphi (x):=[\phi (x_1)\cdots \phi (x_{n_k}){]}^{\top },x\in {\mathbb{R}}^{n_k},$$

$${\left[\begin{array}{c}x\\ 1\end{array}\right]}^{\top }P\left[\begin{array}{c}x\\ 1\end{array}\right]\ge 0\text{for all }x\in \mathcal{X}.$$

$$\mathcal{X}\subseteq \bigcap _{P\in {\mathcal{P}}_{\mathcal{X}}}\left\{x\in {\mathbb{R}}^{n_x}:{\left[\begin{array}{c}x\\ 1\end{array}\right]}^{\top }P\left[\begin{array}{c}x\\ 1\end{array}\right]\ge 0\right\}.$$

$${\mathcal{P}}_{\mathcal{X}}=\left\{P\mid P=\left[\begin{array}{cc}-2\Gamma & \Gamma (\underline{x}+\bar{x})\\ (\underline{x}+\bar{x}{)}^{\top }\Gamma & -2{\underline{x}}^{\top }\Gamma \bar{x}\end{array}\right]\right\},$$

$${\mathcal{P}}_{\mathcal{X}}=\left\{P\mid P=\left[\begin{array}{cc}{H}^{\top }\Gamma H& -H^{\top }\Gamma h\\ -h^{\top }\Gamma H& h^{\top }\Gamma h\end{array}\right]\right\},$$

$${\left[\begin{array}{cc}A& x_c\\ 0& 1\end{array}\right]}^{\top }P\left[\begin{array}{cc}A& x_c\\ 0& 1\end{array}\right]+\left[\begin{array}{cc}2\Gamma & -\Gamma 1_m\\ -1_m^{\top }\Gamma & 0\end{array}\right]⪰0,$$

$${\mathcal{P}}_{\mathcal{X}}=\left\{P\mid P=\mu \left[\begin{array}{cc}-A^{\top }A& -A^{\top }b\\ -b^{\top }A& 1-b^{\top }b\end{array}\right],\mu \ge 0\right\}.$$

$${\mathcal{S}}_x=\bigcap _{i=1}^m\left\{x\in {\mathbb{R}}^{n_x}\mid {\left[\begin{array}{c}x\\ f(x)\\ 1\end{array}\right]}^{\top }{S}_i\left[\begin{array}{c}x\\ f(x)\\ 1\end{array}\right]\le 0\right\},$$

$$S_i=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& c_i\\ 0& c_i^{\top }& -2d_i\end{array}\right]i=1,\cdots ,m.$$

$${\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]}^{\top }Q\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]\ge 0\text{for all }x\in \mathcal{X},$$

$$\mathcal{G}(\varphi )\subseteq \bigcap _{Q\in {\mathcal{Q}}_{\varphi }}\left\{(x,y)\in {\mathbb{R}}^{2n}:{\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]}^{\top }Q\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]\ge 0\right\}.$$

$$(\phi (x)-\alpha x)(\phi (x)-\beta x)\le 0.$$

$$(\varphi (x)-K_{1}x{)}^{\top }(\varphi (x)-K_{2}x)\le 0,$$

$${\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]}^{\top }\left[\begin{array}{ccc}-K_1^{\top }{K}_2-K_2^{\top }{K}_1& K_1^{\top }+K_2^{\top }& 0\\ K_1+K_2& -2I_n& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]\ge 0.$$

$$(\varphi (x)-\varphi (x^{\star })-\alpha (x-x^{\star }){)}^{\top }(\varphi (x)-\varphi (x^{\star })-\beta (x-x^{\star }))\le 0.$$

$$\alpha \le \frac{\varphi (x)-\varphi (x^{\star })}{x-x^{\star }}\le \beta \forall x,x^{\star }\in \mathbb{R}.$$

$$(\phi (x_i)-\phi (x_i^{\star })-\alpha (x_i-x_i^{\star }))(\phi (x_i)-\phi (x_i^{\star })-\beta (x_i-x_i^{\star }))\le 0.$$

$$(\phi (x_i)-\phi (x_j)-\alpha (x_i-x_j))(\phi (x_i)-\phi (x_j)-\beta (x_i-x_j))\le 0,$$

$${\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]}^{\top }\left[\begin{array}{ccc}-2\alpha \beta T& (\alpha +\beta )T& 0\\ (\alpha +\beta )T& -2T& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]\ge 0,$$

$$T=\sum _{1\le i<j\le n}{\lambda }_{ij}(e_i-e_j)(e_i-e_j{)}^{\top },{\lambda }_{ij}\ge 0,$$

$${\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]}^{\top }\left[\begin{array}{ccc}0& 0& 0\\ 0& -2D& D(\underline{\varphi }+\bar{\varphi })\\ 0& (\underline{\varphi }+\bar{\varphi }{)}^{\top }D& -2{\underline{\varphi }}^{\top }D\bar{\varphi }\end{array}\right]\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]\ge 0,$$

$$y_i=\max (\alpha x_i,\beta y_i)⟺\{\begin{array}{l}(y_i-\alpha x_i)(y_i-\beta x_i)=0\\ \beta x_i\le y_i\\ \alpha x_i\le y_i.\end{array}$$

$$(y_j-y_i-\alpha (x_j-x_i))(y_j-y_i-\beta (x_j-x_i))\le 0.$$

$$\sum _{i=1}^n\{{\lambda }_i(y_i-\alpha x_i)(y_i-\beta x_i)-{\nu }_i(y_i-\beta x_i)-{\eta }_i(y_i-\alpha x_i)\}+$$

$$\sum _{i\ne j}{\lambda }_{ij}(y_j-y_i-\alpha (x_j-x_i))(y_j-y_i-\beta (x_j-x_i))\le 0,$$

$${\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]}^{\top }\left[\begin{array}{ccc}{Q}_{11}& Q_{12}& Q_{13}\\ Q_{12}^{\top }& Q_{22}& Q_{23}\\ Q_{13}^{\top }& Q_{23}^{\top }& Q_{33}\end{array}\right]\left[\begin{array}{c}x\\ \varphi (x)\\ 1\end{array}\right]\ge 0,$$

Full text

I Introduction

Neural networks have become increasingly effective at many difficult machine-learning tasks. However, the nonlinear and large-scale nature of neural networks makes them hard to analyze and, therefore, they are mostly used as black-box models without formal guarantees. In particular, neural networks are highly vulnerable to attacks, or more generally, uncertainty in their input. In the context of image classification, for example, neural networks can be easily deluded into changing their classification labels by slightly perturbing the input image [1]. Indeed, it has been shown that even imperceptible perturbations in the input of the state-of-the-art neural networks cause natural images to be misclassified with high probability [2]. Input perturbations can be either of an adversarial nature [3], or they could merely occur due to compression, resizing, and cropping [4]. These drawbacks limit the adoption of neural networks in safety-critical applications such as self-driving vehicles [5], aircraft collision avoidance procedures [6], speech recognition, and recognition of voice commands; see [7] for a survey.

Motivated by the serious consequences of the fragility of neural networks to input uncertainties or adversarial attacks, there has been an increasing effort in developing tools to measure or improve the robustness of neural networks. Many results focus on specific adversarial attacks and attempt to harden the network by, for example, crafting hard-to-classify examples [8, 9, 10, 11]. Although these methods are scalable and work well in practice, they still suffer from false negatives. Safety-critical applications require provable robustness against any bounded variations in the input data. As a result, many tools have recently been used, adapted, or developed for this purpose, such as mixed-integer linear programming [12, 13, 14, 15], convex relaxations and duality theory [16, 17, 18], Satisfiability Modulo Theory (SMT) [19], dynamical systems [20, 21], Abstract Interpretation [22, 23], interval-based methods [24, 25, 26, 27, 28]. All these works aim at bounding the worst-case value of a performance measure when their input is perturbed within a specified range.

Our contribution. In this paper, we develop a novel framework based on semidefinite programming (SDP) for safety verification and robustness analysis of neural networks against norm-bounded perturbations in their input. Our main idea is to abstract nonlinear activation functions of neural networks by the constraints they impose on the pre- and post- activation values. In particular, we describe various properties of activation functions using Quadratic Constraints (QCs), such as bounded slope, bounded values, monotonicity, and repetition across layers. Using this abstraction, any property (e.g., safety or robustness) that we can guarantee for the abstracted network will automatically be satisfied by the original network as well. The quadratic form of these constraints allows us to formulate the verification problem as an SDP feasibility problem. Our main tool for developing the SDP is the $\mathcal{S}$-procedure from robust control [29], which allows us to reason about multiple QCs. Our framework has the following notable features:

• •

We use various forms of QCs to abstract any type of activation function.

• •

Our method can capture the cross-coupling between neurons across different layers, thereby reducing conservatism. This feature, which hinges on the assumption that the same activation function is used throughout the entire network (repetition across layers), becomes particularly effective for deep networks.

• •

We can control the trade-off between computational complexity and conservatism by systematically including or excluding different types of QCs.

In this paper, we focus on the neural network verification problem (formally stated in $\S$II-A) but the proposed framework (input-output characterization of neural networks via quadratic constraints) can be adapted to other problems such as sensitivity analysis of neural networks to input perturbations, output reachable set estimation, probabilistic verification, bounding the Lipschitz constant of neural networks, and closed-loop stability analysis.

I-A Related Work

The performance of certification algorithms for neural networks can be measured along three axes. The first axis is the tightness of the certification bounds; the second axis is the computational complexity, and, the third axis is applicability across various models (e.g. different activation functions). These axes conflict. For instance, the conservatism of the verification algorithm is typically at odds with the computational complexity. The relative advantage of any of these algorithms is application-specific. For example, reachability analysis and safety verification applications call for less conservative algorithms, whereas in robust training, computationally fast algorithms are desirable [24, 16].

On the one hand, formal verification techniques such as Satisfiability Modulo (SMT) solvers [30, 31, 32], or integer programming approaches [14, 15] rely on combinatorial optimization to provide tight certification bounds for piece-wise linear networks, whose complexity scales exponentially with the size of the network in the worst-case. A notable work to improve scalability is [15], where the authors do exact verification of piecewise-linear networks using mixed-integer programming with an order of magnitude reduction in computational cost via tight formulations for non-linearities and careful preprocessing.

On the other hand, certification algorithms based on continuous optimization are more scalable but less accurate. A notable work in this category is [16], in which the authors propose a linear-programming (LP) relaxation of piece-wise linear networks and provide upper bounds on the worst-case loss using weak duality. The main advantage of this work is that the proposed algorithm solely relies on forward- and back-propagation operations on a modified network, and thus is easily integrable into existing learning algorithms. In [33], the authors propose an SDP relaxation of one-layer sigmoid-based neural networks based on bounding the worst-case loss with a first-order Taylor expansion. The closest work to the present work is [34], in which the authors propose a semidefinite relaxation (SDR) for certifying robustness of piece-wise linear multi-layer neural networks. This relaxation is based on the so-called “lifting”, where the original problem is embedded in a much larger space. This SDR approach provides tighter bounds than those of [16] but is less scalable. Finally, compared to the SDR method of [34], our SDP framework yield tighter bounds, especially for deep networks, and is not limited to $\mathrm{ReLU}$ networks. Parts of this work, specialized to probabilistic verification, have appeared in the conference paper [35].

The rest of the paper is organized as follows. In $\S$II we formulate the safety verification problem and present the assumptions. In $\S$III, we abstract the problem with Quadratic Constraints (QCs). In $\S$IV we state our main results. In $\S$V, we discuss further utilities of our framework beyond safety verification. In $\S$VI we provide numerical experiments to evaluate the performance of our method and compare it with competing approaches. Finally, in $\S$VII we draw conclusions.

I-B Notation and Preliminaries

We denote the set of real numbers by $\mathbb{R}$, the set of nonnegative real numbers by $\mathbb{R}_{+}$, the set of real $n$-dimensional vectors by $\mathbb{R}^{n}$, the set of $m\times n$-dimensional matrices by $\mathbb{R}^{m\times n}$, the $m$-dimensional vector of all ones by $\mathrm{1}_{m}$, the $m\times n$-dimensional zero matrix by $0_{m\times n}$, and the $n$-dimensional identity matrix by $I_{n}$. We denote by $\mathbb{S}^{n}$, $\mathbb{S}_{+}^{n}$, and $\mathbb{S}_{++}^{n}$ the sets of $n$-by-$n$ symmetric, positive semidefinite, and positive definite matrices, respectively. The $p$-norm ($p\geq 1$) is displayed by $\|\cdot\|_{p}\colon\mathbb{R}^{n}\to\mathbb{R}_{+}$. For $A\in\mathbb{R}^{m\times n}$, the inequality $A\geq 0$ is element-wise. For $A\in\mathbb{S}^{n}$, the inequality $A\succeq 0$ means $A$ is positive semidefinite. For sets $\mathcal{I}$ and $\mathcal{J}$, we denote their Cartesian product by $\mathcal{I}\times\mathcal{J}$. The indicator function of a set $\mathcal{X}$ is defined as $\mathbf{1}_{\mathcal{X}}(x)=1$ if $x\in\mathcal{X}$, and $\mathbf{1}_{\mathcal{X}}(x)=0$ otherwise. For two matrices $A,B$ of the same dimension, we denote their Hadamard product by $A\circ B$. A function $g\colon\mathbb{R}^{n}\to\mathbb{R}$ is $\alpha$-convex ($0\leq\alpha<\infty$) if $g-(\alpha/2)\|\cdot\|_{2}^{2}$ is convex. Furthermore, $g$ is $\beta$-smooth ($0\leq\beta<\infty$) if it is differentiable and $(\beta/2)\|\cdot\|_{2}^{2}-g$ is convex. Finally, if $g$ is $\alpha$-convex and $\beta$-smooth, then

[TABLE]

for all $x,y\in\mathbb{R}^{n}$ [36, Theorem 2.1.12].

II Safety Verification and Robustness Analysis of Neural Networks

II-A Problem Statement

Consider the nonlinear vector-valued function $f\colon\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{f}}$ described by a multi-layer feed-forward neural network. Given a set $\mathcal{X}\subset\mathbb{R}^{n_{x}}$ of possible inputs (e.g., adversarial examples), the neural network maps $\mathcal{X}$ to an output set $\mathcal{Y}$ given by

[TABLE]

The desirable properties that we would like to verify can often be represented by a safety specification set $\mathcal{S}_{y}$ in the output space of the neural network. In this context, the network is safe if the output set lies within the safe region, i.e., if the inclusion $f(\mathcal{X})\subseteq\mathcal{S}_{y}$ holds. Alternatively, we can define $\mathcal{S}_{x}:=f^{-1}(\mathcal{S}_{y})$ as the inverse image of $\mathcal{S}_{y}$ under $f$. Then, safety corresponds to the inclusion $\mathcal{X}\subseteq\mathcal{S}_{x}$.

Checking the condition $\mathcal{Y}\subseteq\mathcal{S}_{y}$, however, requires an exact computation of the nonconvex set $\mathcal{Y}$, which is very difficult. Instead, our interest is in finding a non-conservative over-approximation ${\mathcal{\tilde{Y}}}$ of $\mathcal{Y}$ and verifying the safety properties by checking the condition $\tilde{\mathcal{Y}}\subseteq\mathcal{S}_{y}$. This approach detects all false negatives but also produces false positives, whose rate depends on the tightness of the over-approximation–see Figure 1. The goal of this paper is to solve this problem for a broad class of input uncertainties and safety specification sets using semidefinite programming.

II-A1 Classification Example

Consider a data (e.g., image) classification problem with $n_{f}$ classes, where a feed-forward neural network $f\colon\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{f}}$ takes as input a data point $x$ and returns an $n_{f}$-dimensional vector of scores (or logits) – one for each class. The classification rule is based on assigning $x$ to the class with the highest score. That is, the class of $x$ is given by $C(x)=\mathrm{argmax}_{1\leq i\leq n_{f}}\ f_{i}(x)$. To evaluate the local robustness of the neural network around a correctly-classified point $x^{\star}$, we consider a set $\mathcal{X}$, containing $x^{\star}$, that represents the set of all possible perturbations of $x^{\star}$. In image classification, a popular choice are perturbations in the $\ell_{\infty}$ norm, i.e., $\mathcal{X}=\{x\colon\|x-x^{\star}\|_{\infty}\leq\epsilon\}$, where $\epsilon$ is the maximum perturbation applied to each pixel. Then the classifier is locally robust at $x^{\star}$ if it assigns all the perturbed inputs to the same class as $x^{\star}$, i.e., if $C(x)=C(x^{\star})$ for all $x\in\mathcal{X}$. For this problem, the safe set is the polytope $\mathcal{S}_{y}=\{y\in\mathbb{R}^{n_{f}}\mid y_{i^{\star}}\geq y_{i}\ \text{for all }i\neq i^{\star}\},$ where $i^{\star}=\mathrm{argmax}_{1\leq i\leq n_{f}}\ f_{i}(x^{\star})$ is the class of $x^{\star}$.

II-B Neural Network Model

For the model of the neural network, we consider an $\ell$-layer feed-forward fully-connected neural network $f\colon\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{f}}$ described by the following recursive equations:

[TABLE]

where $x^{0}\!=\!x\in\mathbb{R}^{n_{0}}(n_{0}\!=\!n_{x})$ is the input to the network and $W^{k}\in\mathbb{R}^{n_{k+1}\times n_{k}},\ b^{k}\in\mathbb{R}^{n_{k+1}}$ are the weight matrix and bias vector of the $(k+1)$-th layer. We denote by $n=\sum_{k=1}^{\ell}n_{k}$ the total number of neurons. The nonlinear activation function $\phi$ ($\mathrm{ReLU}$111Rectified Linear Unit., sigmoid, tanh, etc.) is applied coordinate-wise to the pre-activation vectors, i.e., it is of the form

[TABLE]

where $\varphi$ is the activation function of each neuron. The output $f(x)$ depends on the specific application we are considering. For example, in image classification with cross-entropy loss, $f(x)$ represents the logit input to the softmax function; or, in feedback control, $x$ is the input to the neural network controller (e.g., tracking error) and $f(x)$ is the control input to the plant.

III Problem Abstraction via Quadratic Constraints

In this section, our goal is to provide an abstraction of the verification problem described in $\S$II-A that can be converted into a semidefinite program. Our main tool is Quadratic Constraints (QCs), which were first developed in the context of robust control [37] for describing nonlinear, time-varying, or uncertain components of a system. We start with the abstraction of sets using QCs.

III-A Input Set

We now provide a particular way of representing the input set $\mathcal{X}$ that will prove useful for developing the SDP.

Definition 1

Let $\mathcal{X}\subset\mathbb{R}^{n_{x}}$ be a nonempty set. Suppose $\mathcal{P}_{\mathcal{X}}$ is the set of all symmetric indefinite matrices $P$ such that

[TABLE]

We then say that $\mathcal{X}$ satisfies the QC defined by $\mathcal{P}_{\mathcal{X}}$.

Note that by definition, $\mathcal{P}_{\mathcal{X}}$ is a convex cone, i.e., if $P_{1},P_{2}\in\mathcal{P}_{\mathcal{X}}$ then $\theta_{1}P_{1}+\theta_{2}P_{2}\in\mathcal{P}_{\mathcal{X}}$ for all nonnegative scalars $\theta_{1},\theta_{2}$. Furthermore, we can write

[TABLE]

In other words, we can over approximate $\mathcal{X}$ by the intersection of a possibly infinite number of sets defined by quadratic inequalities. We will see in $\S$IV that the matrix $P\in\mathcal{P}_{\mathcal{X}}$ appears as a decision variable in the SDP. In this way, we can optimize the over-approximation of $\mathcal{X}$ to minimize the conservatism of the specific verification problem we want to solve.

Proposition 1

(QC for hyper-rectangle) The hyper-rectangle $\mathcal{X}=\{x$ $\in\mathbb{R}^{n_{x}}\mid\underline{x}\leq x\leq\bar{x}\}$ satisfies the QC defined by

[TABLE]

where $\Gamma\in\mathbb{R}^{n_{x}\times n_{x}}$ is diagonal and nonnegative. For this set, (5) holds with equality.

Proof 1

See Appendix A-A.

Our particular focus in this paper is on perturbations in the $\ell_{\infty}$ norm, $\mathcal{X}=\{x\mid\|x-x^{\star}\|_{\infty}\leq\epsilon\}$, which are a particular class of hyper-rectangles with $\underline{x}=x^{\star}-\epsilon\mathrm{1}$ and $\bar{x}=x^{\star}+\epsilon\mathrm{1}$. We can adapt the result of Proposition 1 to other sets such as polytopes, zonotopes, and ellipsoids, as outlined below. The derivation of the corresponding QCs can be found in Appendix A-B.

III-A1 Polytopes

Let $\mathcal{X}=\{x\in\mathbb{R}^{n_{x}}\mid Hx\leq h\}$ be a polytope, where $H\in\mathbb{R}^{m\times n_{x}},h\in\mathbb{R}^{m}$. Then $\mathcal{X}$ satisfies the QC defined by

[TABLE]

where $\Gamma\in\mathbb{S}^{m},\Gamma\geq 0,\Gamma_{ii}=0$. Furthermore, if the set $\{x\in\mathbb{R}^{n_{x}}\mid Hx\geq h\}$ is empty, then (5) holds with equality.

III-A2 Zonotopes

A zonotope is an affine transformation of the unit cube, $\mathcal{X}=\{x\in\mathbb{R}^{n_{x}}\mid x=x_{c}+A\lambda,\quad\lambda\in[0,1]^{m}\}$, where $A\in\mathbb{R}^{n_{x}\times m}$ and $x_{c}\in\mathbb{R}^{n_{x}}$. Then any $P\in\mathcal{P}_{\mathcal{X}}$ satisfies

[TABLE]

for some diagonal and nonnegative $\Gamma\in\mathbb{R}^{m\times m}$.

III-A3 Ellipsoids

Suppose the input set $\mathcal{X}$ is an ellipsoid defined by $\mathcal{X}=\{x\in\mathbb{R}^{n_{x}}\mid\|Ax+b\|_{2}\leq 1\}$, where $A\in\mathbb{S}^{n_{x}}$ and $b\in\mathbb{R}^{n_{x}}$. Then $\mathcal{X}$ satisfies the QC defined by

[TABLE]

III-B Safety Specification Set

As mentioned in the introduction, the safe set can be characterized either in the output space ($\mathcal{S}_{y}$) or in the input space ($\mathcal{S}_{x}$). In this paper, we consider the latter. Specifically, we assume $\mathcal{S}_{x}$ can be represented (or inner approximated) by the intersection of finitely many quadratic inequalities:

[TABLE]

where the $S_{i}\in\mathbb{S}^{n_{x}+n_{f}+1}$ are given. In particular, this characterization includes ellipsoids and polytopes in the output space. For instance, for an output safety specification set described by the polytope $\mathcal{S}_{y}=\cap_{i=1}^{m}\left\{y\in\mathbb{R}^{n_{f}}\mid c_{i}^{\top}y-d_{i}\leq 0\right\}$, the $S_{i}$’s are given by

[TABLE]

III-C Abstraction of Nonlinearities by Quadratic Constraints

One of the main difficulties in the analysis of neural networks is the composition of nonlinear activation functions. To simplify the analysis, instead of analyzing the network directly, our main idea is to remove the nonlinear activation functions from the network but retain the constraints they impose on the pre- and post-activation signals. Using this abstraction, any property (e.g., safety or robustness) that we can guarantee for the “constrained” network will automatically be satisfied by the original network as well. In the following, we show how we can encode various properties of activation functions (e.g., monotonicity, bounded slope, and bounded values) using QCs. We first provide a formal definition below.

Definition 2** **(QC for functions)

Let $\phi\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ and suppose $\mathcal{Q}_{\phi}\subset\mathbb{S}^{2n+1}$ is the set of all symmetric indefinite matrices $Q$ such that

[TABLE]

where $\mathcal{X}\subseteq\mathbb{R}^{n}$ is a nonempty set. Then we say $\phi$ satisfies the QC defined by $\mathcal{Q}_{\phi}$ on $\mathcal{X}$.

We remark that our definition of a QC slightly differs from the one used in robust control [37], by including a constant in the vector surrounding the matrix $Q$, which allows us to incorporate affine constraints (e.g., bounded nonlinearities). In view of Definition 1, we can interpret (11) as a QC satisfied by the graph of $\phi$, $\mathcal{G}(\phi):=\{(x,y)\mid y=\phi(x),\ x\in\mathcal{X}\}\subset\mathbb{R}^{2n}$, i.e., $Q_{\phi}=\mathcal{P}_{\mathcal{G}(\phi)}$. Therefore, we can write

[TABLE]

In other words, we over-approximate the graph of $\phi$ by a quadratically constrained set.

The derivation of quadratic constraints is function specific but there are certain rules and heuristics that can be used for all of them which we describe below.

III-C1 Sector-Bounded Nonlinearities

Consider the nonlinear function $\varphi\colon\mathbb{R}\to\mathbb{R}$ with $\varphi(0)=0$. We say that $\varphi$ is sector-bounded in the sector $[\alpha,\beta]$ ($\alpha\leq\beta<\infty$) if the following condition holds for all $x\in\mathbb{R}$,222For the case where $\alpha=-\infty$ or $\beta=+\infty$, we define the sector bound inequality as $x(\varphi(x)-\beta x)\leq 0$ and $x(\alpha x-\varphi(x))\leq 0$, respectively.

[TABLE]

Geometrically, this inequality means that the function $y=\varphi(x)$ lies in the sector formed by the lines $y=\alpha x$ and $y=\beta x$ (see Figure 2). As an example, the $\mathrm{ReLU}$ function belongs to the sector $[0,1]$ and in fact, lies on its boundary.

For the vector case, let $K_{1},K_{2}\in\mathbb{R}^{n\times n}$ be two matrices such that $K_{2}-K_{1}$ is symmetric positive semidefinite. We say that $\phi\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ is sector-bounded in the sector $[K_{1},K_{2}]$ if the following condition holds for all $x\in\mathbb{R}^{n}$ [38],

[TABLE]

or, equivalently,

[TABLE]

The sector condition does not impose any restriction on the slope of the function. This motivates a more accurate description of nonlinearities with bounded slope [39].

III-C2 Slope-Restricted Nonlinearities

A nonlinear function $\phi\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ is slope-restricted in the sector $[\alpha,\beta]$ ($\alpha\leq\beta<\infty$), if for any $x,x^{\star}\in\mathbb{R}^{n}$,

[TABLE]

For the one-dimensional case $(n=1)$, (14) states that the chord connecting any two points on the curve of $\phi$ has a slope that is at least $\alpha$ and at most $\beta$:

[TABLE]

Note that a slope-restricted nonlinearity with $\phi(0)=0$ is also sector bounded. Furthermore, if $\phi$ is slope-restricted in $[\alpha,\beta]$, then the function $x\mapsto\phi(x+x^{\star})-\phi(x^{\star})$ belongs to the sector $[\alpha I_{n},\beta I_{n}]$ for any $x^{\star}$. Finally, the gradient of an $\alpha$-convex and $\beta$-smooth function is slope-restricted in $[\alpha,\beta]$.

To connect the results of the previous subsection to activation functions in neural networks, we recall the following result from convex analysis [36].

Lemma 1** **(gradient of convex functions)

Consider a function $g\colon\mathbb{R}^{n}\to\mathbb{R}$ that is $\alpha$-convex and $\beta$-smooth. Then the gradient function $\nabla g\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$ is slope-restricted in the sector $[\alpha,\beta]$.

Notably, all commonly-used activation functions for deep neural networks are gradients of convex functions. Therefore, they belong to the class of slope-restricted nonlinearities, according to Lemma 1. We have the following result.

Proposition 2

The following statements hold true.

• (a)

The $\mathrm{ReLU}$ function $\varphi(x)=\max(0,x),\ x\in\mathbb{R}$ is slope-restricted and sector-bounded in $[0,1]$.

• (b)

The sigmoid function, $\varphi(x)=\frac{1}{1+e^{-x}},\ x\in\mathbb{R}$ is slope-restricted in $[0,1]$.

• (c)

The tanh function, $\varphi(x)=\tanh(x),\ x\in\mathbb{R}$ is slope-restricted and sector-bounded in $[0,1]$.

• (d)

The leaky $\mathrm{ReLU}$ function, $\varphi(x)=\max(ax,x),\ x\in\mathbb{R}$ with $a>0$ is slope-restricted and sector-bounded in $[\min(a,1),\max(a,1)]$.

• (e)

The exponential linear function (ELU), $\varphi(x)=\max(x,a(e^{x}-1)),\ x\in\mathbb{R}$ with $a>0$ is slope-restricted and sector-bounded in $[0,1]$.

• (f)

The softmax function, $\phi(x)=[\frac{e^{x_{1}}}{\sum_{i=1}^{d}e^{x_{i}}},$ $\cdots,$ $\frac{e^{x_{n}}}{\sum_{i=1}^{d}e^{x_{i}}}]^{\top}$, $x\in\mathbb{R}^{n}$ is slope-restricted in $[0,1]$.

In the context of neural networks, our interest is in repeated nonlinearities of the form $\phi(x)=[\varphi(x_{1})\ \cdots\ \varphi(x_{n})]^{\top}.$ Furthermore, the activation values might be bounded from below or above (e.g., the $\mathrm{ReLU}$ function which outputs a nonnegative value). The sector bound and slope restricted inequalities can become too conservative as they do not capture these properties. In the following, we discuss QCs for these properties.

III-C3 Repeated Nonlinearities

Suppose $\varphi\colon\mathbb{R}\to\mathbb{R}$ is slope-restricted in $[\alpha,\beta]$ ($\alpha\leq\beta$) and let $\phi(x)=[\varphi(x_{1})$ $\cdots\varphi(x_{n})]^{\top}$ be a vector-valued function constructed by component-wise repetition of $\varphi$. It is not hard to verify that $\phi$ is also slope-restricted in $[\alpha,\beta]$. Indeed, by summing the slope-restriction conditions

[TABLE]

over $i=1,\cdots,n$, we obtain (14). However, this representation simply ignores the fact that all the nonlinearities that compose $\phi$ are the same. By taking advantage of this structure, we can refine the quadratic constraint that describes $\phi$. To be specific, for an input-output pair $(x,\phi(x)),\ x\in\mathbb{R}^{n}$, we can write the inequality

[TABLE]

for all distinct $i,j=1,\cdots,n,\ i\neq j$. This particular QC can considerably reduce conservatism, especially for deep networks, as it reasons about the coupling between the neurons throughout the entire network. By making an analogy to dynamical systems, we can interpret the neural network as a time-varying discrete-time dynamical system where the same nonlinearity is repeated for all “time” indexes $k$ (the layer number). Then the QC in (15) couples all the possible neurons. In the following lemma, we characterize QCs for repeated nonlinearities.

Lemma 2

(QC for repeated nonlinearities)* Suppose $\varphi:\mathbb{R}\to\mathbb{R}$ is slope-restricted in the sector $[\alpha,\beta]$. Then the vector-valued function $\phi(x)=[\varphi(x_{1})\ \cdots\varphi(x_{n})]^{\top}$ satisfies the QC*

[TABLE]

for all $x\in\mathbb{R}^{n}$, where

[TABLE]

and $e_{i}\in\mathbb{R}^{n}$ is the $i$-th unit vector.

Proof 2

By a conic combination of $\binom{n}{2}$ quadratic constraints of the form (15), we obtain (16). See Appendix A-C for a detailed proof.

There are several results in the literature about repeated nonlinearities. For instance, in [40, 41], the authors derive QCs for repeated and odd nonlinearities (e.g. tanh function).

III-C4 Bounded Nonlinearities

Finally, suppose the nonlinear function values are bounded, i.e., $\underline{\phi}\leq\phi(x)\leq\bar{\phi}$ for all $x\in\mathbb{R}^{n}$. Using Proposition 1, $\phi(x)$ satisfies the quadratic constraint

[TABLE]

for all $x$, where $D\in\mathbb{R}^{n\times n}$ is diagonal and nonnegative. We can write a similar inequality when the pre- activation values are known to be bounded. More generally, if the graph of $\phi$ is known to satisfy $\mathcal{G}(\phi)\subseteq\mathcal{G}$, then any quadratic constraint for $\mathcal{G}$ is also a valid quadratic constraint for $\phi$.

We observe that the inequalities (15)-(18) are all quadratic in $(x,\phi(x),1)$, and therefore can be encapsulated into QCs of the form (11). As we show in $\S$IV, the matrix $Q\in\mathcal{Q}_{\phi}$ that abstracts the nonlinearity $\phi$ appears as a decision variable in the SDP.

Although the above rules can be used to guide the search for valid QCs for activation functions, a less conservative description of activation functions requires a case-by-case treatment to further exploit the structure of the nonlinearity. In the next subsection, we elaborate on QCs for $\mathrm{ReLU}$ activation functions.

III-D Quadratic Constraints for $\mathrm{ReLU}$ Activation Function

The $\mathrm{ReLU}$ function precisely lies on the boundary of the sector $[0,1]$. This observation can be used to refine the QC description of $\mathrm{ReLU}$. Specifically, let $y=\max(\alpha x,\beta x),\ x\in\mathbb{R}^{n}$ be the concatenation of $n$ $\mathrm{ReLU}$ activation functions333For $\mathrm{ReLU}$, we have $\alpha=0$ and $\beta=1$.. Then each individual activation function can be described by the following constraints [34]:

[TABLE]

The first constraint is the boundary of the sector $[\alpha,\beta]$ and the other constraints simply prune these boundaries to recover the $\mathrm{ReLU}$ function. Furthermore, for any two distinct indices $i\neq j$, we can write the constraint (15):

[TABLE]

By adding a weighted combination of all these constraints (non-negative weights for inequalities), we find that the function $y=\max(\alpha x,\beta x)$ satisfies

[TABLE]

for all $x\in\mathbb{R}^{n}$. In the following lemma, we provide a full QC characterization of the $\mathrm{ReLU}$ function.

Lemma 3** **(Global QC for $\mathrm{ReLU}$ function)

The function $\phi(x)=\max(\alpha x,\beta x)$ satisfies the QC

[TABLE]

for all $x\in\mathbb{R}^{n}$, where

[TABLE]

$\nu,\eta\in\mathbb{R}_{+}^{n}$*, and $T$ is given by (17). *

Proof 3

See Appendix A-D.

III-D1 Tightening Relaxations

The QC of Lemma 3 holds globally for the whole space $\mathbb{R}^{n}$. When restricted to a local region $\mathcal{X}$, these QCs can be tightened. Specifically, suppose $y=\max(x,0)$ and define $\mathcal{I}^{+}$, $\mathcal{I}^{-}$, and $\mathcal{I}^{\pm}$ as the set of activations that are known to be always active, always inactive, or unknown for all $x\in\mathcal{X}\subseteq\mathbb{R}^{n}$, i.e.,

[TABLE]

Then the function $y_{i}=\max(\alpha x_{i},\beta x_{i})$ belongs to the sector $[\alpha,\alpha],\ [\alpha,\beta]$ and $[\beta,\beta]$ for inactive, unknown, and active neurons, respectively. Furthermore, since the constraint $y_{i}\geq\beta x_{i}$ holds with equality for active neurons, we can write $\nu_{i}\in\mathbb{R}\text{ if }i\in\mathcal{I}^{+},\ \nu_{i}\geq 0\text{ otherwise}.$ Similarly, the constraint $y_{i}\geq\alpha x_{i}$ holds with equality for inactive neurons. Therefore, we can write $\eta_{i}\in\mathbb{R}\text{ if }i\in\mathcal{I}^{-},\ \eta_{i}\geq 0\text{ otherwise}$. Finally, the chord connecting the input-output pairs of always-active or always-inactive neurons has slope of $\alpha$ or $\beta$. Equivalently, for any $(i,j)\in(\mathcal{I}^{+}\times\mathcal{I}^{+})\cup(\mathcal{I}^{-}\times\mathcal{I}^{-})$, we can write

[TABLE]

Therefore, in (III-D), $\lambda_{ij}\in\mathbb{R}$ for $(i,j)\in(\mathcal{I}^{+}\times\mathcal{I}^{+})\cup(\mathcal{I}^{-}\times\mathcal{I}^{-})$ and $\lambda_{ij}\geq 0$ otherwise. The above additional degrees of freedom on the multipliers can tighten the relaxation incurred in (III-D). In the following Lemma, we summarize the above observations.

Lemma 4

(Local QC for $\mathrm{ReLU}$ function) Let $\phi(x)=\max(\alpha x,\beta x),\ x\in\mathcal{X}\subset\mathbb{R}^{n}$ and define $\mathcal{I}^{+},\mathcal{I}^{-}$ as in (23). Then $\phi$ satisfies the QC

[TABLE]

for all $x\in\mathcal{X}$, where

[TABLE]

with $T=\sum_{1\leq i<j\leq n}\lambda_{ij}(e_{i}-e_{j})(e_{i}-e_{j})^{\top}$ and

[TABLE]

Proof 4

See Appendix A-D.

We do not know a priori which neurons are always active or always inactive. However, we can partially find them by computationally cheap presolve steps. Specifically, if $x$ is known to satisfy $\underline{x}\leq x\leq\bar{x}$ (bounds on the pre-activation values), then we have $\mathcal{I}^{+}=\{i\mid\underline{x}_{i}\geq 0\}$, $\mathcal{I}^{-}=\{i\mid\bar{x}_{i}<0\}$, and $\mathcal{I}^{\pm}=\{i\mid\bar{x}_{i}\underline{x}_{i}\leq 0\}$. These element-wise bounds can be found by, for example, interval bound propagation [42, 43] or the LP approach of [16]. Indeed, tighter bounds result in a less conservative description of the $\mathrm{ReLU}$ function outlined in Lemma 4.

III-E Other Activation Functions

Deriving non-conservative QCs for other activation functions (other than $\mathrm{ReLU}$) is more complicated as they are not on the boundary of any sector. However, by bounding these functions at multiple points by sector bounds of the form (12), we can obtain a substantially better over-approximation. In Figure 3, we illustrate this idea for the $\tanh$ function.

A secondary approach is to use the element-wise bounds on the inputs to the activation functions to use a tighter sector bound condition in (12). For instance, suppose $x\in[\underline{x},\bar{x}]\subseteq\mathbb{R}$. Then the function $\varphi(x)=\tanh(x)$ satisfies the sector condition in (12), where $\alpha$ and $\beta$ are given by

[TABLE]

More generally, suppose the graph of $\varphi\colon[\underline{x},\bar{x}]\to\mathbb{R}$ is known to satisfy $\mathcal{G}(\varphi)\subseteq\mathcal{G}\subset\mathbb{R}^{2}$. Then any QC satisfied by $\mathcal{G}$ is also a valid QC for $\varphi$. We can use this property to build local quadratic constraints for general activation functions provided that we can overapproximate their graph locally. This idea is illustrated in Figure 3 for the case of $\tanh$ function.

IV Neural Network Verification Via Semidefinite Programming

In the previous section, we developed an abstraction of sets and nonlinearities using quadratic constraints. In this section, we use this abstraction to develop an LMI feasibility problem that can assert whether $f(\mathcal{X})\subseteq\mathcal{S}_{y}$ (or $\mathcal{X}\subseteq\mathcal{S}_{x}=f^{-1}(\mathcal{S}_{y})$). The crux of our idea in the development of the LMI is the $\mathcal{S}$-procedure [29], a technique to reason about multiple quadratic constraints, and is frequently used in robust control and optimization [44, 45].

IV-A Single-layer Neural Networks

For the sake of simplicity in the exposition, we start with the analysis of one-layer neural networks and then extend the results to the multi-layer case in $\S$IV-B. We further assume that the safe set $\mathcal{S}_{x}$ in (10) is specified by a single quadratic form, i.e., $m=1$. We state our main result in the following theorem.

Theorem 1** **(SDP for one layer)

Consider a one-layer neural network $f\colon\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{f}}$ described by the equation

[TABLE]

Suppose $x\in\mathcal{X}\subset\mathbb{R}^{n_{x}}$, where $\mathcal{X}$ satisfies the QC defined by $\mathcal{P}_{\mathcal{X}}$, i.e., for any $P\in\mathcal{P}_{\mathcal{X}}$,

[TABLE]

Let $\mathcal{Z}=\{z\mid z=W^{0}x+b^{0},\ x\in\mathcal{X}\}$ and suppose $\phi$ satisfies the QC defined by $\mathcal{Q}_{\phi}$ on $\mathcal{Z}$, i.e., for any $Q\in\mathcal{Q}_{\phi}$,

[TABLE]

*Consider the following matrix inequality, *

[TABLE]

where

[TABLE]

and $S\in\mathbb{S}^{n_{x}+n_{f}+1}$ is a given symmetric matrix. If (28) is feasible for some $P\in\mathcal{P}_{\mathcal{X}},\ Q\in\mathcal{Q}_{\phi}$, then

[TABLE]

Proof 5

See Appendix A-F.

Theorem 1 states that if the matrix inequality (28) is feasible for some $(P,Q)\in\mathcal{P}_{\mathcal{X}}\times\mathcal{Q}_{\phi}$, then we can certify that the network $\mathcal{X}\subseteq\mathcal{S}_{x}$ or $f(\mathcal{X})\subseteq\mathcal{S}_{y}$. Since $\mathcal{P}_{\mathcal{X}}$ and $\mathcal{Q}_{\phi}$ are both convex, (28) is a linear matrix inequality (LMI) feasibility problem and, hence, can be efficiently solved via interior-point method solvers for convex optimization.

Remark 1** **(End-to-end QC for neural network)

It follows from the proof of Theorem 1 that, in view of Definition 2, the neural network in (25) satisfies the QC defined by $(\mathcal{X},\mathcal{Q}_{f})$, where

[TABLE]

In other words, for any $Q_{f}\in\mathcal{Q}_{f}$ we have

[TABLE]

IV-B Multi-layer Neural Networks

We now turn to multi-layer neural networks. Assuming that all the activation functions are the same across the layers (repetition across layers), we can concatenate all the pre- and post-activation signals together and form a more compact representation. To see this, we first introduce ${\mathbf{x}}=[{x^{0}}^{\top}\cdots{x^{\ell}}^{\top}]^{\top}\in\mathbb{R}^{n_{0}+n}$, where $\ell\geq 1$ is the number of hidden layers. We further define the entry selector matrices $\mathbf{E}^{k}\in\mathbb{R}^{n_{k}\times(n_{0}+n)}$ such that $x^{k}=\mathbf{E}^{k}\mathbf{x}$ for $k=0,\cdots,\ell$. Then, we can write (2) compactly as

[TABLE]

In the following result, we develop the multi-layer counterpart of Theorem 1 for the multi-layer neural network in (31).

Theorem 2** **(SDP for multiple layers)

Consider the multi-layer neural network described by (31). Suppose $\mathcal{X}\subset\mathbb{R}^{n_{x}}$ satisfies the QC defined by $\mathcal{P}_{\mathcal{X}}$. Define $\mathcal{Z}=\{\mathbf{A}\mathbf{x}+\mathbf{b}\mid x\in\mathcal{X}\}$ and suppose $\phi$ satisfies the QC defined by $\mathcal{Q}_{\phi}$ on $\mathcal{Z}$. Consider the following LMI.

[TABLE]

where

[TABLE]

and $S\in\mathbb{S}^{n_{x}+n_{f}+1}$ is a given symmetric matrix. If (32) is feasible for some $(P,Q)\in\mathcal{P}_{\mathcal{X}}\times\mathcal{Q}_{\phi}$, then

[TABLE]

Proof 6

See Appendix A-G.

Remark 2

For the case that the safe set is characterized by more than one quadratic inequality, i.e., when $m>1$ in (10), then $\mathcal{X}\subseteq\mathcal{S}_{x}$ if the following LMIs,

[TABLE]

hold for some $P_{i}\in\mathcal{P}_{\mathcal{X}}$ and $Q_{i}\in\mathcal{Q}_{\phi}$.

V Optimization Over the Abstracted Network

In the previous section, we developed an LMI feasibility problem as a sufficient to verify the safety of the neural network. We can incorporate this LMI as a constraint of an optimization problem to solve problems beyond safety verification. Specifically, we can define the following SDP,

[TABLE]

where $g(P,Q,S)$ is a convex function of $P,Q,S$, and $\mathcal{S}$ is a convex subset of $\mathbb{S}^{n_{x}+n_{f}+1}$. In the following, we allude to some utilities of the SDP (36), which we call $\mathrm{DeepSDP}$.

V-A Reachable Set Estimation

In Theorem 1, we developed a feasibility problem to assert whether $\mathcal{X}\subseteq\mathcal{S}_{x}$, or equivalently, $f(\mathcal{X})\subseteq\mathcal{S}_{y}$. By parameterizing $\mathcal{S}_{x}$, we can find the best over approximation of $f(\mathcal{X})$ by solving (36). Suppose $\mathcal{S}_{x}$ is described by $\mathcal{S}_{x}=\{x\mid c^{\top}f(x)-d\leq 0\}$ with a given $c\in\mathbb{R}^{n_{f}}$ and $d\in\mathbb{R}$. By defining

[TABLE]

the feasibility of (32) for some $(P,Q)\in\mathcal{P}_{\mathcal{X}}\times\mathcal{Q}_{\phi}$ implies $c^{\top}f(x)\leq d\quad\text{for all }x\in\mathcal{X}$. In other words, $d$ is a certified upper bound on the optimal value of the optimization problem

[TABLE]

Now if we treat $d\in\mathbb{R}$ as a decision variable, we can minimize this bound by solving (36) with $g(P,Q,S)=d$. This is particularly useful for over approximating the reachable set $f(\mathcal{X})$ by a polyhedron of the form ${\mathcal{S}_{y}}=\cap_{i}\left\{y\in\mathbb{R}^{n_{f}}\mid c_{i}^{\top}y-d_{i}\leq 0\right\}$, where $c_{i}$ are given and the goal is to find the smallest value of $d_{i}$, for each $i$, such that $f(\mathcal{X})\subseteq{\mathcal{S}_{y}}$.

By reparameterizing $S$ in (37) we can also compute the best ellipsoidal over-approximation of $f(\mathcal{X})$. Specifically, define

[TABLE]

Then the inclusion $f(\mathcal{X})\subseteq\mathcal{S}_{y}=f(\mathcal{S}_{x})$ implies that $f(\mathcal{X})$ is enclosed by the ellipsoid $\mathcal{S}_{y}=\{y\in\mathbb{R}^{n_{f}}\mid\|A_{y}y+b_{y}\|_{2}\leq 1\}$. Therefore, finding the minimum-volume ellipsoid enclosing $f(\mathcal{X})$ amounts to the optimization problem

[TABLE]

Note that this problem is not convex in $(A_{y},b_{y})$ due to the non-affine dependence of $S$ on these variables. However, by using Schur complements, we can formulate an equivalent convex program. We skip the details for the sake of space and refer the reader to [35].

V-B Closed-Loop Reachability Analysis

By modifying the matrix $S$ in (37), we can use a similar approach as presented in $\S\ref{subsec: Certified Upper Bounds}$ to over approximate the reachable sets of closed-loop systems involving neural networks. Specifically, consider a discrete-time Linear Time-Invariant (LTI) system driven by a neural network controller,

[TABLE]

Given a set of current states $\mathcal{X}$, the one-step forward reachable set is $\mathcal{X}^{+}=f_{cl}(\mathcal{X}).$ Suppose $\mathcal{S}_{x}$ in (10) is defined by $\mathcal{S}_{x}=\{x\in\mathbb{R}^{n_{x}}\mid c^{\top}f_{cl}(x)\leq d\}$, where

[TABLE]

According to Theorem 2, the feasibility of the LMI (32) for some $(P,Q)\!\in\!\mathcal{P}_{\mathcal{X}}\times\mathcal{Q}_{\phi}$ would allow us to conclude $\mathcal{X}\subseteq\mathcal{S}_{x}$, or equivalently, $c^{\top}f_{cl}(x)\leq d$ for all $x\in\mathcal{X}$. By repeating this for different pairs $(c_{i},d_{i})\in\mathbb{R}^{n_{f}}\times\mathbb{R},\ i=1,\cdots,m$, we can overapproximate the one-step reachable set $f_{cl}(\mathcal{X})$ by the polyhedron $\mathcal{P}=\left\{x^{+}\in\mathbb{R}^{n_{x}}\mid c_{i}^{\top}x^{+}-d_{i}\leq 0\ i=1,\cdots,m\right\}$. Similarly, we can also overapproximate the closed-loop reachable sets by ellipsoids. In $\S\ref{subsection: Verification of Approximate Model Predictive Control}$ we use this approach to verify a model predictive controller approximated by a neural network.

VI Discussion and Numerical Experiments

In this section, we discuss the numerical aspects of our approach. For solving the SDP, we used MOSEK [46] with CVX [47] on a 5-core personal computer with 8GB of RAM. For all experiments, we used $\mathrm{ReLU}$ activation functions and did Interval Bound Propagation as a presolve step to determine the element-wise bounds on the activation functions444All code, data, and experiments for this paper are available at https://github.com/mahyarfazlyab/DeepSDP. . We start with the computational complexity of the proposed SDP.

VI-A Computational Complexity

VI-A1 Input Set

The number of decision variables for the input set depends on the set representation. The quadratically constrained set that over-approximates hyperectangles is indexed by $n_{x}$ decision variables, where $n_{x}$ is the input dimension (see Proposition 1). Note that for hyper-rectangles, we can include additional quadratic constraints. Indeed, any $x$ satisfying $\underline{x}\leq x\leq\bar{x}$ satisfies $2n_{x}^{2}-n_{x}$ quadratic constraints of the form $(x_{i}-\underline{x}_{i})(\bar{x}_{j}-x_{i})\geq 0$, $(x_{i}-\underline{x}_{i})(x_{j}-\underline{x}_{j})\geq 0\ i\neq j$, $(x_{i}-\bar{x}_{i})(x_{j}-\bar{x}_{j})\geq 0\ i\neq j$. However, one can precisely characterize a hyper-rectangle with only $n_{x}$ of these quadratic constraints, namely, $(x_{i}-\underline{x}_{i})(\bar{x}_{i}-x_{i})\geq 0$. Our numerical computations reveal that adding the remaining QCs would not tighten the relaxation.

For polytopes, the maximum number of decision variables is $\binom{m}{2}$, where $m$ is the number of half-spaces defining the polytope. However, we can use some heuristics to remove quadratic constraints that are not “tight”. For instance, for the polytope $\mathcal{X}=\{x\mid Hx\leq h\}$, we can write $\binom{m}{2}$ sector bounds of the form $(H_{i}^{\top}x-h_{i})(H_{j}^{\top}x-h_{j})\geq 0.$ Now if the intersection of these hyperplanes belongs to $\mathcal{X}$, then the sector would be tight (see Figure 5). We can verify this by checking the feasibility of

[TABLE]

Finally, for the case of ellipsoids, we only have one decision variable, the parameter $\mu$ in (9).

VI-A2 Activation Functions

For a network with $n$ hidden neurons, if we use all possible quadratic constraints, the number of decision variables will be $\mathcal{O}(n+n^{2})$. If we ignore repeated nonlinearities, we will arrive at $\mathcal{O}(n)$ decision variables. In our numerical experiments, we did not observe any additional conservatism after removing repeated nonlinearities across the neurons of the same layer. However, accounting for repeated nonlinearities was sometimes very effective for the case of multiple layers.

VI-A3 Safety Specification Set

The number of decision variables for the safety specification set depends on how we would like to bound the output set. For instance, for finding a single hyperplane, we have only one decision variable. For the case of ellipsoids, there will be $\mathcal{O}(n_{f}^{2})$ decision variables.

VI-B Synthetic Examples

VI-B1 Number of Hidden Layers

As the first experiment, we consider finding over-approximations of the reachable set of a neural network with a varying number of layers, for a given input set. Specifically, we consider randomly-generated neural networks with $n_{x}=2$ inputs, $n_{f}=2$ outputs, and $\ell=\{1,2,3,4\}$ hidden layers, each having $n_{k}=100$ neurons per layer. For the input set, we consider $\ell_{\infty}$ balls with center $x^{\star}=(1,1)$ and radius $\epsilon=0.1$. We use $\mathrm{DeepSDP}$ to find over-approximations of $f(\mathcal{X})$ in the form of polytopes (see $\S$V-A). In Figure 4, we compare the output set $f(\mathcal{X})$ (using exhaustive search over $\mathcal{X}$) with two over-approximations: the red polytope is obtained by solving $\mathrm{DeepSDP}$. The dashed black polytope is obtained by the semidefinite relaxation ($\mathrm{SDR}$) approach of [34]. We observe that the bounds obtained by $\mathrm{DeepSDP}$ are relatively tighter, especially for deeper networks. In Appendix A-H, we provide more visualizations.

VI-B2 Repeated Nonlinearities

As the second experiment, we study the effect of including repeated nonlinearities on the tightness of the bounds. Specifically, we bound the output of a randomly-generated neural network with $n_{x}=2$ inputs, $n_{f}=2$ output, and $n_{k}=10$ neurons per layer by a polytope with 6 facets. For the input set we consider $\ell_{\infty}$ ball with center $x^{\star}=(1,1)$ and radius $\epsilon=0.1$. In Figure 6, we plot the output set, and it over-approximation by $\mathrm{DeepSDP}$ before and after including repeated nonlinearities. We observe that by including repeated nonlinearities, the bounds become tighter, especially for deep networks.

VI-B3 Comparison with Other Methods

As the third experiment, we consider the following optimization problem,

[TABLE]

To evaluate the tightness of our bounds, we compare $\mathrm{DeepSDP}$ with the $\mathrm{MILP}$ formulation of [13], the semidefinite relaxation ($\mathrm{SDR}$) of [34], and the LP relaxation of [16]. For the problem data, we generated random instances of neural networks with $n_{x}=10$ inputs, $n_{f}=1$ output and $\ell\in\{1,\cdots,5\}$ hidden layers; for each layer size, we generated $100$ random neural networks with their weights and biases chosen independently from the normal distribution $\mathcal{N}(0,1/\sqrt{n_{x}})$. For the input set, we consider $x^{\star}=\mathrm{1}_{n_{x}}$ and $\epsilon=0.2$. In Table I, we report the comparisons of bounds and running times. The MILP formulation finds the global solution but the running time grows quickly as the number of neurons increases. Compared to $\mathrm{SDR}$, the bounds of $\mathrm{DeepSDP}$ are relatively tighter, especially for deeper networks. Finally, the LP relaxation bounds are considerably looser but the running time is negligible. In Figure 7, we plot the histograms of the normalized gap between the optimal value $f^{\star}$ (obtained by MILP) and the upper bound $f^{\mathrm{SDP}}$ for layer sizes $\ell=1,2$.

VI-C Verification of Approximate Model Predictive Control

Consider an LTI system

[TABLE]

where $x_{k}\in\mathbb{R}^{n_{x}}$ is the state at time $k$, $u_{k}\in\mathbb{R}^{n_{u}}$ is the control input, and $A,B$ are matrices of appropriate size. The state and control input are subject to the box constraints $\mathcal{X}=\{x\mid\underline{x}\leq x\leq\bar{x}\}$ and $\mathcal{U}=\{u\mid\underline{u}\leq u\leq\bar{u}\}$.

Suppose the control policy is parameterized by a multi-layer fully-connected feed-forward network $f$ that is trained off-line to approximate a model predictive control (MPC) law $\mu^{\star}(x)$. The motivation is to reduce the computational burden of solving an optimization problem online to determine the MPC control action. The trained neural network, however, does not necessarily satisfy the specifications of the MPC control law such as state and control constraint satisfaction. To ensure input constraint satisfaction, we project the neural network output onto $\mathcal{U}$, resulting in the closed-loop system,

[TABLE]

Note that for input box constraints, $\mathcal{U}=\{u\mid\underline{u}\leq u\leq\bar{u}\}$, we can embed the projection operator as two additional layers with a specific choice of weights and biases. Indeed, for an $\ell$-layer $f$, we can describe $f_{p}(x)=\mathrm{Proj}_{\mathcal{U}}(f(x_{k}))$ via the $(\ell+2)$-layer $\mathrm{ReLU}$ network,

[TABLE]

To validate state constraint satisfaction, we must ensure that there is a set of initial states $\mathcal{E}\subseteq\mathcal{X}$ whose trajectories would always satisfy the state constraints. One such set is a positive invariant set. By definition, a set $\mathcal{E}$ is positively invariant with respect to $f_{cl}$, if and only if $x_{0}\in\mathcal{E}$ implies $x_{k}\in\mathcal{E}$ for all $k\geq 1$. Equivalently, $\mathcal{E}$ is positively invariant if $f_{cl}(\mathcal{E})\subseteq\mathcal{E}$. We now show that how we can compute a positive invariant set for (43) using semidefinite programming.

To find a positive invariant set for the closed-loop system, we consider the candidate set $\mathcal{E}=\{x\mid\|x\|_{\infty}\leq\epsilon\}$. We first over approximate the one-step reachable set $f_{cl}(\mathcal{E})$ by the polytope $\mathcal{P}=\{x\mid Hx\leq h\}$, $H\in\mathbb{R}^{m\times n_{x}},h\in\mathbb{R}^{m}$ (see $\S\ref{subsection:Closed-Loop Reachability Analysis}$). To do this, we form the following $m$ SDPs

[TABLE]

where

[TABLE]

With this choice of $S_{i}$, it is not difficult to show that the feasibility of the LMIs in (45) implies $f_{cl}(\mathcal{E})\subseteq\mathcal{P}$, and therefore, (45) finds the smallest $\mathcal{P}$ that encloses $f_{cl}(\mathcal{E})$. Then, $\mathcal{E}$ is positively invariant if $\mathcal{P}\subseteq\mathcal{E}$.

For the numerical experiment, we first consider a 2D system

[TABLE]

subject to the state and input constraints $x_{t}\in\mathcal{X}=\{x\mid\|x\|_{\infty}\leq 5\}$ and $u\in\mathcal{U}=\{u\mid\|u\|_{\infty}\leq 1\}$. We are interested in stabilizing the system by solving the finite horizon problem

[TABLE]

and choosing the control law as $\mu_{MPC}(x)=u_{0}^{\star}$. For generating the training data, we compute $\mu_{MPC}(x)$ at 6284 uniformly chosen random points from the control invariant set. We then train a neural network with two inputs, one output, and two hidden layers with 32 and 16 neurons, respectively using the mean-squared loss. In Figure 9, we plot the explicit MPC control law as well as its approximation by the neural network.

In Figure 8, we plot the largest invariant set $\mathcal{E}$ that we could find, which is $\mathcal{E}=\{x\mid\|x\|_{\infty}\leq 0.65\}$. In this figure, we also plot the output reachable sets for the first four time steps, starting from the initial set $\mathcal{E}$, as well as their over-approximations by $\mathrm{DeepSDP}$.

VII Conclusions

We proposed a semidefinite programming framework for robustness analysis and safety verification of feed-forward fully-connected neural networks with general activation functions. Our main idea is to abstract the nonlinear activation functions by quadratic constraints that are known to be satisfied by all possible input-output instances of the activation functions. We then showed that we can analyze the abstracted network via semidefinite programming. We conclude this paper with several future directions.

First, a notable advantage of the proposed SDP compared to other convex relaxations is the relative tightness of the bounds. In particular, coupling all pairs of neurons in the network (repeated nonlinearities) can considerably reduce conservatism. However, coupling all neurons is not feasible for even medium-sized networks as the number of decision variables would scale quadratically with the number of neurons. Nevertheless, our numerical experiments show that most of these pair-wise couplings of neurons are redundant and do not tighten the bounds. It would be interesting to develop a method that can decide a priori that coupling which pairs of neurons would tighten the relaxation. Second, one of the drawbacks of SDPs is their limited scalability in general. Exploiting the structure of the problem (e.g. sparsity patterns induced by the network strucrure) to reduce the computational complexity would be an important future direction. Third, we have only considered fully-connected networks in this paper. It would be interesting to extend the results to other architectures. Finally, incorporating the proposed framework in training neural networks with desired robustness properties would be another important future direction.

Appendix A Appendix

A-A Proof of Proposition 1

The inequality $\underline{x}\leq x\leq\bar{x}$ is equivalent to $n_{x}$ quadratic inequalities of the form $(x_{i}-\underline{x}_{i})(\bar{x}_{i}-x_{i})\geq 0\quad i=1,\cdots,n_{x}$. Multiplying both sides of with $\Gamma_{i}\geq 0$, summing over $i=1,\cdots,n_{x}$, and denoting $\Gamma=\operatorname{diag}(\gamma_{1},\cdots,\gamma_{n_{x}})$ yields the claimed inequality. $\square$

A-B QCs for Polytopes, Zonotopes, and Ellipsoids

A-B1 Polytopes

For every vector $x$ satisfying $Hx\leq h$, we have $(H_{i}^{\top}x-h_{i})(H_{j}^{\top}x-h_{j})\geq 0,\ i\neq j$, where $H_{i}^{\top}$ is the $i$-th row of $H$. These inequalities imply

[TABLE]

where $\Gamma_{ij}=\Gamma_{ji}\geq 0,\ i\neq j,\ \Gamma_{ii}=0$. The preceding inequality is equivalent to (7). Now suppose the set $\{x\mid Hx\geq h\}$ is empty. Then

[TABLE]

To show this set equality define $\mathcal{X}_{Q}$ as the set on the right-hand side. We have $\mathcal{X}\subset\mathcal{X}_{Q}$. To show $\mathcal{X}_{Q}\subset\mathcal{X}$, suppose $x\in\mathcal{X}_{Q}$, implying that either $H_{i}^{\top}x-h_{i}\leq 0$ for all $i$ or $H_{i}^{\top}x-h_{i}\geq 0$ for all $i$. But the latter cannot happen since the set $\{x\mid Hx\geq h\}$ is empty. Therefore, we have $H_{i}^{\top}x-h_{i}\leq 0\text{ for all }i$.

A-B2 Zonotopes

By multiplying both sides of (8) by $[\lambda^{\top}\ 1]$ and $[\lambda^{\top}\ 1]^{\top}$, respectively, and noting that $x=x_{c}+A\lambda$ we obtain

[TABLE]

where the right inequality follows from the fact that $\lambda\in[0,1]^{m}$, hence satisfying the QC of Proposition 1.

A-B3 Ellipsoids

Any $x\in\mathcal{X}$ satisfies $\mu(1-(Ax+b)^{\top}(Ax+b))\geq 0$ for $\mu\geq 0$. The latter inequality is equivalent to (9). $\square$

A-C Proof of Lemma 2

For any distinct pairs $(x_{i},\varphi(x_{i}))$ and $(x_{j},\varphi(x_{j}))$, $1\leq i<j\leq n$, we can write the slope restriction inequality in (14) as

[TABLE]

By multiplying both sides by $\lambda_{ij}\geq 0$, we obtain

[TABLE]

where $E_{ij}=(e_{i}-e_{j})(e_{i}-e_{j})^{\top}$ and $e_{i}\in\mathbb{R}^{n}$ is the $i$-th unit vector in $\mathbb{R}^{n}$. Summing over all $1\leq i<j\leq n$ will yield the desired result.

A-D Proof of Lemma 3

Consider the equivalence in (19) for the $i$-th coordinate of $y=\max(\alpha x,\beta x),\ x\in\mathbb{R}^{n}$:

[TABLE]

Multiplying these constraints by $\lambda_{i}\in\mathbb{R}$, $\nu_{i}\in\mathbb{R}_{+}$, and $\eta_{i}\in\mathbb{R}_{+}$, respectively, and adding them together, we obtain

[TABLE]

Substituting $x_{i}=e_{i}^{\top}x$ and $y_{i}=e_{i}^{\top}y$, where $e_{i}$ is the $i$-th unit vector in $\mathbb{R}^{n}$, and rearranging terms, we get

[TABLE]

where

[TABLE]

Furthermore, since $y_{i}=\max(\alpha x_{i},\beta x_{i})$ is slope-restricted in $[\alpha,\beta]$, by Lemma 2 we can write

[TABLE]

Summing (48) over all $i=1,\cdots,n$ and adding the result to (49) would yield (22). $\square$

A-E Proof of Lemma 4

Consider the relation $y=\max(\alpha x,\beta x)$. For active neurons, $i\in\mathcal{I}^{+}$, we can write

[TABLE]

Similarly, for inactive neurons, $i\in\mathcal{I}^{-}$,we can write

[TABLE]

Finally, for unknown neurons, $i\in\mathcal{I}^{\pm}$, we can write

[TABLE]

A weighted combination of the above constraints yields

[TABLE]

where $\alpha_{i}=\alpha+(\beta-\alpha)\mathbf{1}_{\mathcal{I}^{+}}(i)$, $\beta_{i}=\beta-(\beta-\alpha)\mathbf{1}_{\mathcal{I}^{-}}(i)$, $\nu_{i}\in\mathbb{R}_{+}\text{ for }i\notin\mathcal{I}^{+}$ and $\eta_{i}\in\mathbb{R}_{+}\text{ for }i\notin\mathcal{I}^{-}$. Furthermore, since $y_{i}=\max(\alpha x_{i},\beta x_{i})$ is slope-restricted on $[\alpha,\beta]$, we can write

[TABLE]

Adding (50) and (51) and rearranging terms would yield the desired inequality. $\square$

A-F Proof of Theorem 1

Consider the identity $x^{1}=\phi(W^{0}x^{0}+b^{0})$. Using the assumption that $\phi$ satisfies the quadaratic constraint defined by $\mathcal{Q}_{\phi}$ on $\mathcal{Z}$, $x^{0},x^{1}$ satisfy the QC

[TABLE]

for any $Q\in\mathcal{Q}_{\phi}$ and all $x^{0}\in\mathcal{X}$. By assumption $\mathcal{X}$ satisfies the QC defined by $\mathcal{P}_{\mathcal{X}}$, implying that for any $P\in\mathcal{P}_{\mathcal{X}}$,

[TABLE]

for all $x^{0}\in\mathcal{X}$. Suppose (28) holds for some $(P,Q)\in\mathcal{P}_{\mathcal{X}}\times\mathcal{Q}_{\phi}$. By left- and right- multiplying both sides of (28) by $[{x^{0}}^{\top}\ {x^{1}}^{\top}\ 1]$ and $[{x^{0}}^{\top}\ {x^{1}}^{\top}\ 1]^{\top}$, respectively, we obtain

[TABLE]

Therefore, the last term on the left-hand side must be nonpositive for all $x^{0}\in\mathcal{X},\ x^{1}=\phi(W^{0}x^{0}+b^{0})$, or, equivalently,

[TABLE]

Using the relations $x^{0}=x$ and $f(x)=W^{1}x^{1}+b^{1}$, the above inequality is the desired inequality in (34). $\square$

A-G Proof of Theorem 2

Recall the definition $\mathcal{Z}=\{\mathbf{A}\mathbf{x}+\mathbf{b}\mid x\in\mathcal{X}\}$. Since $\phi$ satisfies the QC defined by $\mathcal{Q}_{\phi}$ on $\mathcal{Z}$, for any $Q\in\mathcal{Q}_{\phi}$, we have

[TABLE]

By assumption $\mathcal{X}$ satisfies the QC defined by $\mathcal{P}_{\mathcal{X}}$. Using the relation $x^{0}=\mathbf{E}^{0}{\mathbf{x}}$, for any $P\in\mathcal{P}_{\mathcal{X}}$ it holds that

[TABLE]

Suppose the LMI in (32) holds for some $(P,Q)\in\mathcal{P}_{\mathcal{X}}\times\mathcal{Q}_{\phi}$. By left- and right- multiplying both sides of (27) by $[{\mathbf{x}}^{\top}\ 1]$ and $[{\mathbf{x}}^{\top}\ 1]^{\top}$, respectively, we obtain

[TABLE]

Therefore, the last quadratic term must be nonpositive for all $x^{0}\in\mathcal{X}$, from where we can write

[TABLE]

Using the relations $x^{0}=\mathbf{E}^{0}{\mathbf{x}}$ and $f(x)=W^{\ell}\mathbf{E}^{\ell}{\mathbf{x}}+b^{\ell}$ from (31), the above inequality can be written as

[TABLE]

$\square$

A-H More Visualizations

In Figure 10, we show the effect of the number of hidden neurons on the quality of approximation for a single-layer network, and in Figure 11, we change the perturbation size.