POLAR: A POLYNOMIAL ARITHMETIC FRAMEWORK FOR VERIFYING NEURAL-NETWORK CONTROLLED SYSTEMS
Year: 2021

Abstract:
We propose POLAR, a polynomial arithmetic framework that leverages polynomial overapproximations with interval remainders for bounded-time reachability analysis of neural-network controlled systems (NNCSs). Compared with existing arithmetic approaches that use standard Taylor models, our framework uses a novel approach to iteratively overapproximate the neuron output ranges layer-by-layer via a combination of Bernstein polynomial interpolation for continuous activation functions and Taylor model arithmetic for the other operations. This approach overcomes the main drawback in the standard Taylor model arithmetic, i.e. its inability to handle functions that cannot be well approximated by Taylor polynomials, and significantly improve the accuracy and efficiency of reachable states computation for NNCSs. To further tighten the overapproximation, our method keeps the Taylor model remainders symbolic under the linear mappings when propagating Taylor models across the neural-network controller. We show that POLAR can be seamlessly integrated with existing Taylor model flowpipe construction techniques, and our approach significantly outperforms the current state-of-the-art techniques on a suite of benchmarks.

heading: 1 INTRODUCTION
text: Neural networks have been increasingly used as the central decision makers in a variety of tasks Mnih et al. (2015); Lillicrap et al. (2016); Pan et al. (2018). However, the use of neural-network controllers also gives rise to new challenges on verifying the correctness of the resulting closed-loop control systems especially in safety-critical settings. In this paper, we consider the reachability verification problem of neural-network controlled systems (NNCSs). The high-level architecture of a simple NNCS is shown in Figure 1 in which the neural network senses the system state, i.e. the value of x, at discrete time steps, and computes the corresponding control values u for updating the system dynamics which is defined by an ordinary differential equation (ODE) over x and u. The boundedtime reachability analysis problem of an NNCS is to compute an (over-approximated) reachable set that contains all the trajectories starting from the initial set for a finite number of control steps. Figure 2 shows an illustration of reachable sets for 4 steps, where the orange regions are the reachable set, and two red arrowed curves are two exemplary trajectories starting from the initial set X0 (blue).
Reachability analysis of general NNCSs is notoriously difficult due to nonlinearity in both the neural-network controller and the plant, further exacerbated by the coupling of the controller and the plant. Since exact reachability of general nonlinear systems is undecidable Alur & Dill (1994), current approaches for reachability analysis of nonlinear dynamical systems aim at computing a tight overapproximation of the reachable sets Dreossi et al. (2016); Lygeros et al. (1999); Yang et al. (2016); Prajna & Jadbabaie (2004); Huang et al. (2017a); Frehse et al. (2011); Chen et al. (2013); Althoff (2015). Earlier works on NNCS verification drew on techniques for computing the output ranges of neural networks Huang et al. (2017b); Katz et al. (2017); Dutta et al. (2018); Wang et al. (2018); Weng et al. (2018); Zhang et al. (2018); Singh et al. (2019) when integrating with the aforementioned reachability analysis. However, they have been shown to be ineffective for NNCS verification due to the lack of consideration on the interaction between the neural-network controller and the plant dynamics Dutta et al. (2019); Huang et al. (2019). In particular, since their primary goal is to bound the output range of the neural network instead of approximating its input-output function,
Figure 1: A typical NNCS model.
they cannot be used to track dependency across the closed-loop system and across multiple time steps in reachability analysis.
Direct end-to-end approximation Dutta et al. (2019); Huang et al. (2019); Fan et al. (2020) and layer-by-layer propagation Ivanov et al. (2019; 2021b;a) are the two main categories of approaches in the NNCS verification literature. The former directly computes a function approximation of the neural network and suffers from efficiency problems, i.e. they cannot handle systems with more than a few dimensions due to the need to sample in the input space. The latter approach tries to exploit the neural network structure and uses Taylor model (TM) arithmetic to more efficiently obtain a function approximation by propagating the TMs layer by layer. A Taylor model (p, I) consists of a polynomial p for point-wise approximation, and an interval remainder I to bound the approximation error. However, due to the limitations of basic TM arithmetic, these approaches cannot handle non-differentiable activation functions and suffer from rapid growth of the interval remainders during propagation which effectively degrades the analysis to an interval analysis.
In this paper, we tackle the challenge of dependency tracking by constructing a function overapproximation, specifically a Taylor model (TM) approximation, of the neural-network controller using a unified polynomial arithmetic framework (POLAR) that enables precise layer-by-layer propagation of TMs for general feed-forward neural networks. POLAR addresses the key challenges of applying basic TM arithmetic through a novel use of univariate Bernstein polynomial interpolation to handle non-differentiable activation functions and symbolic remainders to taper the growth of the interval remainders. For the former, we leverage an efficient sampling-based analysis to provide a sound overapproximation. The latter avoids the so-called wrapping effect Jaulin et al. (2001) that can lead to rapid growth of the interval remainders in linear mappings.
In summary, our work makes the following novel contributions.
• We propose a novel polynomial arithmetic framework for bounded-time reachability analysis of NNCSs, keeping track of state-wise dependency across the closed-loop system. Compared to existing propagation-based approaches, our framework has the advantage of being able to handle NN controllers with general and heterogeneous activation functions.
• We propose neuron-wise Bernstein polynomial interpolation and show that it can be seamlessly integrated with Taylor model approximations. In addition, we present the first application of symbolic remainders to tightening the overapproximation of neural network behaviors.
• We conduct a comprehensive comparison of our approach with state-of-the-art techniques, including an evaluation on the full set of benchmarks published in these works, showing the overwhelming advantage of our proposed approach.

heading: 2 PRELIMINARIES
text: A Neural-Network Controlled System (NNCS) is a continuous plant governed by a neural network controller. The plant dynamics is defined by an ODE of the form ẋ = f(x,u) wherein the state variables and control inputs are denoted by the vectors x and u respectively. We assume that the function f is at least locally Lipschitz continuous such that its solution w.r.t. an initial state and constant control inputs is unique Meiss (2007). We denote the input-output mapping of the neural network controller as κ, and the controller is triggered every δc time which is also called the control
stepsize, then a system execution (trajectory) is produced in the following way. Starting from an initial state x(0), the controller senses the system state at the beginning of every control step t=jδc for j=0, 1, . . ., and update the control inputs to vj=κ(x(jδc)). The system’s dynamics in that control step is governed by the ODE ẋ=f(x,vj). Given an initial state set X0 ⊂ Rn, all executions from a state in it can be formally defined by a flowmap function ϕN : X0 × R≥0 → Rn, such that the system state at any time t ≥ 0 from any initial state x0 ∈ X0 is ϕN (x0, t). We call a state x′ ∈ Rn reachable if there exists x0 ∈ X0 and t ≥ 0 such that x′ = ϕN (x0, t). The reachability problem on NNCS is to decide whether a state is reachable in a given NNCS, and it is undecidable due to their higher expressiveness than the two-counter machines on which the reachability problem is undecidable Alur & Dill (1994). Many formal verification problems can be reduced to the reachability problem. For example, the safety verification problem can be reduced to the reachability problem to an unsafe state. In the paper, we focus on computing the reachable set for an NNCS over a bounded number K of control steps. Since flowmap ϕN often does not have a closed form due to the nonlinear ODEs, we seek to compute state-wise overapproximations for it over time segments, that is, in each control step [(j − 1)δc, jδc] for j = 1, . . . ,K, the reachable set is overapproximated by a group of flowpipes F1(x0, τ), . . . ,FN (x0, τ) over the N uniformly subdivided time segments of the time interval, such that Fi(x0, τ) is a state-wise overapproximation of ϕN (x0, (j − 1)δc + (i− 1)δ + τ) for τ ∈ [0, δc/N ], i.e., Fj(x0, τ) contains the exact reachable state from any initial state x0 in the i-th time segment of the j-th control step. Here, τ is the local time variable which is independent in each flowpipe. A high-level flowpipe construction algorithm is presented as follows, in which X̂0 = X0 and δ = δc/N which is called the flowpipe step or time step.
1: for j = 1 to K do 2: Computing the an overapproximation Ûj−1 for the control input range κ(X̂j−1); 3: Computing the flowpipes F1(x0, τ), . . . ,FN (x0, τ) for the continuous dynamics ẋ = f(x,u), u̇ = 0 from the initial set x(0) ∈ X̂j−1, u(0) ∈ Ûj−1; 4: R ← R∪ {(F1(x0, τ), . . . ,FN (x0, τ)}; 5: X̂j ← FN (z, δ); 6: end for
Notice that x(0) denotes the local initial set for the ODE used in the current control step, while the variables x0 in a flowpipe are the symbolic representation of an initial state in X0. Intuitively, a flowpipe overapproximates not only the reachable set in a time step, but also the dependency from an initial state to its reachable state at a particular time.
In the paper, we use Taylor models to represent the flowpipes. Taylor models are originally proposed to compute higher-order overapproximations for the ranges of continuous functions (see Berz & Makino (1998)). They can be viewed as a higher-order extension of intervals Moore et al. (2009) which are sets of real numbers between lower and upper real bounds, e.g., the interval [a, b] wherein a ≤ b represents the set of {x | a ≤ x ≤ b}. A Taylor model (TM) is a pair (p, I) wherein p is a polynomial of degree k over a finite group of variables x1, . . . , xn ranging in an interval domain D ⊂ Rn, and I is the remainder interval. The range of a TM is the Minkowski sum of the range of its polynomial and the remainder interval. Thereby we sometimes intuitively denote a TM (p, I) by p+ I in the paper. TMs are often used as overapproximations for smooth functions. Given a smooth function f(x) with x ∈ D, its order k TM overapproximation, or simply TM, can be obtained as (pf (x), If ) such that pf is the order k Taylor approximation of f w.r.t. a point in D, and If is a remainder interval which ensures that ∀x ∈ D.(f(x) ∈ pf (x) + If ). Since pf is also over the same variables x, the overestimation in a TM can be estimated only on the width of If , i,e., if the remainder is zero, then pf defines the same mapping as f . At this point, a TM serves more like an overapproximate function than just a range overapproximation such as intervals or polyhedra.
TM arithmetic. TMs are closed under operations such as addition, multiplication, and integration (see Makino & Berz (2003)). Given functions f, g that are overapproximated by TMs (pf , If ) and (pg, Ig) respectively. A TM for f + g can be computed as (pf + pg, If + Ig), and an order k TM for f · g can be computed as ( pf · pg − rk , If · B(pg) + B(pf ) · Ig + If ·Ig + B(rk) ) wherein B(p) denotes an interval enclosure of the range of p, and the truncated part rk consists of the terms in pf · pg of degrees > k. Similar to reals and intervals, TMs can also be organized as vectors and matrices to overapproximate the functions whose ranges are multidimensional. For a vector of TMs ((p1, I1), . . . , (pn, In))
T such that p1, . . . , pn are over the same variables, we collectively denote it by (p, I) such that p is the polynomial vector (p1, . . . , pn)T and I is interval vector (I1, . . . , In)T .
As an example, given the TMs (1− 0.5x2, [−0.1, 0.1]) and (x+ 0.1x4, [−0.2, 0.2]) over the domain x ∈ [−1, 1]. The order 4 TM for the sum is (1 + x − 0.5x2 + 0.1x4, [−0.3, 0.3]), and the order 4 TM for the product is (x− x3 + x4, [−0.38, 0.38]). Although the technique of Taylor model flowpipe construction can be used to generate state-wise overapproximate flowpipes for ODEs (see Berz & Makino (1998); Chen (2015)), it is still a challenge that how to compute and represent the set Ûj−1 by a TM which overapproximates the dependency from an initial state inX0 to the actual control input rangeUj−1 = {κ(ϕN (x0, (j−1)δc)) |x0 ∈ X0} at the beginning of that control step. We illustrate this dependency in Figure 3. Our approach aims at solving this problem.

heading: 3 POLYNOMIAL ARITHMETIC FOR NNCS VERIFICATION
text: In this section, we show how to compute a TM representation for the set Û to overapproximate the dependency from an initial state in X0 to the actual control input range U in terms of the neural network controller κ.
Obviously, the key challenges of providing a TM overapproximation (pr(x0), Ir) in terms of κ lie on how to properly construct a polynomial pr to well capture the state-wise dependency of κ, and how to effectively estimate a remainder Ir to provide a sound and tight error bound. While the basic TM arithmetic that works in a layer-by-layer manner for function approximation of a neural network shows a great potential, demonstrated in early works Ivanov et al. (2021b;a), it still suffers from two main challenges in terms of the construction of pr and estimation of Ir respectively.
1. Construction of pr. Taylor approximation requires the approximated function to be differentiable, thus excluding the use of common operations such as ReLU and pooling, and
2. Estimation of Ir. Interval remainders grows rapidly during layer-by-layer propagation of the TMs, effectively degrading it to an interval analysis.
Our approach, POLAR, focuses exactly on these two challenges by a novel use of univariate Bernstein polynomial interpolation (BP) to handle non-differentiable activation functions and symbolic remainders to taper the growth of the interval remainders.
Meanwhile, it is worthy noting that Bernstein polynomial interpolation can also help improve the remainder estimation, as it is always more accurate than Taylor expansion at the same order for a function. Thus, we use univariate Bernstein polynomial interpolation for all the activation functions, not limited to non-differentiable ReLU.
Algorithm 1 Polynomial arithmetic for overapproximating neural network outputs Input: Input TM (p(x0), I) with x0 ∈ X0, the M + 1 matrices W1, . . . ,WM+1 of the weights on
the incoming edges of the hidden and the output layers, the M + 1 vectors B1, . . . , BM+1 of the neurons’ bias in the hidden and the output layers, the M + 1 activation functions σ1, . . . , σM+1 of hidden and output layers. Output: a TM (pr(x0), Ir) that contains the set κ((p(x0), I)). 1: (pr, Ir)← (p, I); 2: for i = 1 to M + 1 do 3: (pt, It) ← Wi · (pr, Ir) +Bi; # Using TM arithmetic 4: Computing an order k Bernstein polynomial pσ for the activation function over the range
defined by (pt, It), i.e., pσ(y) ≈ σi(y) with y ∈ (pt, It); 5: Evaluating a safe remainder Iσ for the polynomial pσ such that ∀y ∈ (pt, It).(σi(y) ∈ pσ(y) + Iσ) holds; 6: (pr, Ir) ← pσ(pt + It) + Iσ; # Using TM arithmetic 7: end for 8: return (pr, Ir).
We use the following running example to demonstrate the effectiveness of POLAR. We show the result of using basic TM arithmetic here, and in each of the following sections, we will demonstrate the improvements from using Bernstein polynomials and symbolic remainders respectively.
Example 1 (Running example with basic TM arithmetic) We demonstrate the computation of a TM output for a small neural network shown in Figure 4. The input is a vector of two TMs over x0 = (z1, z2)
T which denotes the initial state in a 2-D state space. The neural network has 2 inputs, 1 output and 2 hidden layers each of which has 2 neurons. The activation functions of all the hidden layers and output layer are sigmoid. In the figure, the TMs 1− 0.5z1 + z2 − 0.3z1z2 + [−0.1, 0.1] and −2 + z2 − 0.1z21 + [−0.2, 0.2] are the inputs, and z1, z2 ∈ [−1, 1]. We hope to obtain an order 2 TM for the output. Note that, with the basic TM arithmetic, we cannot obtain the output TM, as the remainders for the first hidden layer already explode to the order of 1e23 due to the coarse Taylor approximations for the neural network.

heading: 3.1 POLYNOMIAL ARITHMETIC INTEGRATED WITH BERNSTEIN POLYNOMIAL
text: Given a TM input for a neural network, our approach combines the use of TM arithmetic and Bernstein polynomial interpolation to produce a state-wise overapproximation set for the output of a neural network. The high-level steps are described in Algorithm 1, in which we assume that the neural network has M hidden layers, the neurons in the same layer have the same type of activation functions, and the output layer’s neurons also have activation functions. We collectively use σ(y) to denote applying the activation function σ to each element of y.
In Algorithm 1, the input TM is propagated layer-wisely. While the linear transformation in terms of the layer weight and bias (line 3) and the TM composition (line 6) can be calculated by TM arithmetic, the key step is how to point-wisely approximate activation function by Bernstein polynomial (line 4) and how to evaluate a sound remainder efficiently (line 5).
Computing the Bernstein polynomial pσ (line 4). Given the computed TM (pt, It) in Line 3, pσ is a vector of univariate Bernstein polynomials Lorentz (2013) each component of which is expressed in a variable yj which denotes the range of the j-th component (dimension) of (pt, It). Then the
j-th component of pσ is computed as pjσ(yj) = ∑k s=0 ( σ( Ȳj−Y j k s+ Y j) ( k s ) (yj−Y j)s(Ȳj−yj)k−s (Ȳj−Y j)k ) which is an approximation of the σ(yj) such that Ȳj and Y j denotes the upper and lower bounds respectively of the range in the j-th dimension in (pt, It), and they can be obtained by interval evaluating the TM.
Evaluating a remainder Iσ for pσ (line 5). The remainder Iσ is a vector of intervals, each component of which is a conservative remainder for the corresponding Bernstein polynomial in pσ such that the approximation error for σ((pt, It)) is contained. To do so, the j-th interval in Iσ is evaluated as
[− j , j ] such that
j = max i=1,··· ,m (∣∣∣∣∣pjσ(Y j − Y jm (i− 12) + Y j)− σj(Y j − Y jm (i− 12) + Y j) ∣∣∣∣∣+Lj · Y j − Y jm ) wherein m is the number of samples which are uniformly selected to estimate the remainder. The soundness of the error bound estimation above has been proven in Huang et al. (2019) for multivariate Bernstein polynomials. Since univariate Bernstein polynomials, which we use in this paper, is a special case of multivariate Bernstein polynomials, our approach is also sound. We also provide a detailed proof in the appendix.
Advantages of Bernstein polynomials. We briefly discuss the advantages of using BP in the approximation of activation functions. (i) Taylor approximation requires the approximated function to be differentiable, however BP approximation only requires the function to be continuous. This fact makes Taylor approximation not applicable to ReLU functions. (ii) BP approximations are essentially polynomial interpolations which are always more accurate approximation forms than Taylor expansion at the same order for the a function (see Phillips (2003)). Thus, we apply univariate BP for general activation functions, rather than just for ReLU (iii) Unlike using multivariate Bernstein polynomials in Huang et al. (2019), our framework only computes univariate Bernstein polynomials such that the complexity is always linear in the approximation order k, and the remainder is much easier to be evaluated. In contrast, the size of a multivariate Bernstein polynomial grows exponentially in k.
Example 2 (Running example with Bernstein polynomial) We can obtain the output TM of Example 1 with BP (see the supplementary document), where the remainder is [−0.132, 0.132].
It is important that the polynomial arithmetic algorithm 1 containing the univariate BP can provide a sound state-wise approximation of a neural network, that leads to the soundness of overall flowpipe construction. Specifically, if the input TM (p, I) is a state-wise overapproximation for the reachable state at the beginning of i-th control step, i.e., p(x0) + I contains the reachable state at the time t = (i− 1)δc from x0 for all x0 ∈ X0, then the returned TM pr(x0) + Ir by Algorithm 1 contains the exact i-th control input in the execution starting from x0 for all x0 ∈ X0, i.e., the output TM (pr, Ir) is a state-wise overapproximation of the control input. Hence, we have the important property given in Theorem 1, and then the setR computed by flowpipe construction (Section 2) is an overapproximation of the system reachable set in K control steps. A detailed proof can be found in the appendix.
Theorem 1 If (p(x0, t), I) is the i-th flowpipe computed in the j-th control step, then for any initial state x0 ∈ X0, the box (p(x0, t), I) is guaranteed to contain the reachable state ϕN (x0, (j − 1)δc + (i− 1)δ + t), i.e., the TM flowpipes computed by our approach are state-wise overapproximations for the reachable sets.

heading: 3.2 SYMBOLIC REMAINDERS IN NEURAL NETWORK ANALYSIS
text: The layer-by-layer propagation across the neural network in Algorithm 1 requires a sequence of consecutive compositions of TMs (line 3 and line 6 for each layer), then the TM remainders may potentially be greatly enlarged by these operations, since they are evaluated by pure interval arithmetic. For example, we want to obtain a TM for the polynomial p(y) = c+Ay+ph(y) with y = (q(x0), I), such that c+Ay is the constant and linear part of p, and ph is the higher-order part. The result can be obtained by evaluating c + A · q(x0) + A · I + ph(q(x0) + I) by TM arithmetic such that the term A · I will be wrapped by its interval evaluation which introduces extra overestimation. During the consecutive compositions of TMs, that overestimation will accumulate, and such phenomenon is called wrapping effect Jaulin et al. (2001). To further reduce the overestimation in our TMs, we propose to adapt the symbolic remainder technique, which was previously proposed in Chen & Sankaranarayanan (2016) for ODEs, to avoid the wrapping effect in the layer-by-layer propagation for NNs under linear mappings. The key idea is that we do not compute the interval of A · I out when propagating the TMs to the next layer, but keeps its linear transformation matrix A. When the set is involved in a later linear mapping B ·A · I , we only composing the matrices and keep the set as (B ·A) · I . Such a technique requires us to use extra lists to keep the linear transformation matrices. The detailed method is given below.
Algorithm 2 TM output computation using symbolic remainders, input and output are the same as those in Algorithm 1
1: (pr, Ir)← (p, I); 2: Setting Q as an empty array which can keep M + 1 matrices; 3: Setting J as an empty array which can keep M + 1 multidimensional intervals; 4: for i = 1 to M + 1 do 5: Computing the composite function qi and the remainder interval Iσi using Bernstein overapproximation and TM arithmetic (similar to Line 3 to 5 in Algorithm 1); 6: Constructing Qi using the coefficients of the linear terms in qi; 7: Computing qRi from qi by eliminating the linear terms; 8: (φi, Ji)← qRi (pr + Ir) + Iσi ; # Using TM arithmetic 9: J← Ji; 10: for j = 1 to i− 1 do 11: Q[j]← Qi · Q[j]; 12: end for 13: Adding Qi to Q as the last element; 14: for j = 2 to i do 15: J← J +Q[j] · J [j − 1]; 16: end for 17: Adding Ji to J as the last element; 18: Computing (pr, Ir) according to (1); # Using TM arithmetic 19: end for 20: return R.
Specifically, we denote the resulting TM (pr, Ir) and the Bernstein overapproximation (pσ, Iσ) in the i-th iteration (layer) in Algorithm 1 by (pr,i, Ir,i) and (pσi , Iσi) respectively. Then (pr,i(x0), Ir,i) for i = 1, . . . ,M + 1 is computed as pr,i(x0) + Ir,i = pσi(Wi · (pr,i−1(x0) + Ir,i−1) + Iσi with (pr,0(x0), Ir,0) = (p(x0), I) which is the input TM. If we use qi(y) to denote the function pσi(Wi ·y), the above expression can be simplified to pr,i(x0)+Ir,i = qi(pr,i−1(x0)+Ir,i−1)+Iσi . By decomposing the linear and the remaining part in qi, the function can be expressed as qi(y) = Qiy + q R i (y) such that Qi is the constant matrix consists of the coefficients of the linear terms in qi, and qRi is the nonlinear and the constant part in qi. Hence, the output TM can be further expressed as
pr,i(x0) + Ir,i = Qi(pr,i−1(x0) + Ir,i−1) + q R i (pr,i−1(x0) + Ir,i−1) + Iσi︸ ︷︷ ︸
φi(x0)+Ji
and then by unfolding the recurrence relation, we have
pr,i(x0) + Ir,i = Qi(pr,i−1(x0) + Ir,i−1) + φi(x0) + Ji
= Qi(Qi−1(pr,i−2(x0) + Ir,i−2) + φi−1(x0) + Ji−1) + φi(x0) + Ji = Qi · · · · ·Q1 · p(x0) +Qi · · · · ·Q1 · I + Φi(x0) + Ji (1)
wherein Φi = φi(x0) + Qi · φi−1(x0) + · · · + Qi · · · · · Q2 · φ1(x0), and Ji = Ji + Qi · Ji−1 + · · ·+Qi · · · · ·Q2 · J1. Hence, the term Qi · · · · ·Q1 · I can be evaluated without wrapping effect, and if we keep an array for Jj , the computation of Ji is also free from wrapping effect. Algorithm 2 shows the improvement of Algorithm 1 in accuracy using the symbolic remainder method such that we additionally use two arrays: Q[j] is Qi · · · · ·Qj and J [j] is Ji for 1 ≤ j ≤ i.
Example 3 (Running example with symbolic remainder) The remainder of the output TM of Example 1 with our symbolic remainder technique is further reduced to [−0.069, 0.069], which is more than 40% improvement compared to the result without symbolic remainder, i.e. [−0.132, 0.132].
Time and space complexity. Although Algorithm 2 always produces a TM with a smaller remainder than Algorithm 1 because of the symbolic treatment of the remainder intervals I and Ji under linear mappings, it requires (1) two extra arrays to keep the matrices QM , QM ·QM−1, . . . , QM · · · · ·Q1, and the remainder JM+1, . . . , J1, (2) two extra inner loops which perform i− 1 and i− 2 iterations in the i-th outer iteration. The size of Qi · · · · ·Qj is determined by the rows in Wi and the columns in Wj , and hence the maximum number of neurons in a layer determines the maximum size of the
matrices in Q. Similarly, the maximum dimension of Ji is also bounded by the maximum number of neurons in a layer. Because of the two inner loops, time complexity of Algorithm 2 is quadratic in M , whereas Algorithm 1 is linear in M .

heading: 4 EXPERIMENTS
text: In this section, we first present an illustrating example of attitude control with 6 state variables and 3 control inputs. We then present a comprehensive comparison to the state-of-the-art tools over the full benchmarks in the related work Huang et al. (2019). Finally, we remark on the observed limitations of our approach. All our experiments were run on a machine with 6-core 2.90 GHz Intel Core i5 and 8GB of RAM.
Illustrating example: attitude control. We consider the attitude control of a rigid body with six states and three inputs as a physically illustrating example Prajna et al. (2004). The complexity of this example lies in the combination of the numbers of the state variables and control inputs. The system dynamics is ω̇1 = 0.25(u0 + ω2ω3), ω̇2 = 0.5(u1 − 3ω1ω3), ω̇3 = u2 + 2ω1ω2, ψ̇1 = 0.5 ( ω2(ψ 2 1 + ψ 2 2 + ψ 2 3 − ψ3) + ω3(ψ21 + ψ22 + ψ2 + ψ23) + ω1(ψ21 + ψ22 + ψ23 + 1) ) , ψ̇2 = 0.5 ( ω1(ψ 2 1 + ψ 2 2 + ψ 2 3 + ψ3) + ω3(ψ 2 1 − ψ1 + ψ22 + ψ23) + ω2(ψ21 + ψ22 + ψ23 + 1) ) ,
ψ̇3 = 0.5 ( ω1(ψ 2 1 + ψ 2 2 − ψ2 + ψ23) + ω2(ψ21 + ψ1 + ψ22 + ψ23) + ω3(ψ21 + ψ22 + ψ23 + 1) ) .
wherein the state x=(ω, ψ) consists of the angular velocity vector in a body-fixed frame ω∈R3, and the Rodrigues parameter vector ψ∈R3. The control torque u∈R3 is updated every 0.1 seconds by a neural network with 3 hidden layers each of which has 64 neurons. The activation of the hidden layers are sigmoid and identity respectively. The initial state set is: ω1∈[−0.45,−0.44], ω2∈[−0.55,−0.54], ω3∈[0.65, 0.66], ψ1∈[−0.75,−0.74], ψ2 ∈ [0.85, 0.86], ψ3 ∈ [−0.65,−0.64]. POLAR computed the TM flowpipes for 30 control steps in 201 seconds. Figure 5 shows the plot of the octagonal enclosures of the flowpipes. We can observe that the flowpipes computed by POLAR are tight w.r.t. the simulated traces. As a comparison, although Verisig 2.0 Ivanov et al. (2021a) can handle this system in theory, its remainder explodes very quickly and the tool crashes after only a few steps.
Comparison to state-of-the-arts. We compare POLAR to the state-of-the-art tools on the full benchmarks in Huang et al. (2019), including Sherlock Dutta et al. (2019) (only for ReLU) , Verisig 2.0 Ivanov et al. (2021a) (only for sigmoid and tanh), NNV Tran et al. (2020), and ReachNN*Fan et al. (2020)1. We refer to Huang et al. (2019); Ivanov et al. (2021a) for the details of the benchmarks. The results are presented in Table 1 where NNV is not included since we were not able to successfully use it to prove any of the benchmarks likely because it is designed for linear systems. Similar results for NNV are also observed in Ivanov et al. (2021a). We can see that POLAR successfully verifies almost all cases within seconds, which is on average 23x and up to 353x faster compared to the tool with the second best efficiency, and achieves the best performance among all the tools (detailed
1The results of ReachNN* are based on GPU acceleration.
comparisons on the tightness of reachable sets can be found in the Appendix). In addition, POLAR scales better with the size of the neural network controller compared to ReachNN* and Verisig 2.0. This demonstrates the potential of handling larger scale systems such as the attitude control example. Note that the only exception is #2 with sigmoid NN controller. We discuss this in detail below.
Limitations. POLAR did not prove the property on test # 2 with sigmoid activation but Verisig 2.0 and ReachNN* did. The main reason is the combined use of Bernstein and Taylor approximations. Although a lower-order Taylor approximation can be derived from a higher-order one by truncating the higher-order terms, it is not the case for a Bernstein approximation since the term coefficients often need to be recomputed if the order is changed. Thus, simply using truncation to simplify a TM may produce a large remainder. Secondly, TM-based techniques can be sensitive to the choice of TM settings (stepsizes, orders, etc.). In our experiments, we use the same TM settings for all the tools.

heading: 5 CONCLUSION
text: In this paper, we propose POLAR, a polynomial arithmetic framework, which integrates TM flowpipe construction, Bernstein overapproximation, and symbolic remainder method to efficiently compute reachable set overapproximations for NNCS. Empirical comparison over a suite of benchmarks show that POLAR performs significantly better than state-of-the-art techniques in terms of both computation efficiency and tightness of reachable set estimation. Future work includes automatic tuning of the hyper-parameters in POLAR to further improve performance and accessibility.
Reproducibility Statement. The implementation code can be found in our attached attached zip package, and we will also release it on GitHub once published. In addition, we also provide proofs of two important propositions, the TM outputs of Example 2 and Example 3, and the additional experimental results that include the visualization of the reachable set computed by different approaches in the Appendix.

heading: A TM OUTPUTS OF THE RUNNING EXAMPLE
text: The result with Bernstein polynomial:
(7.920880969687617e− 01) + (−7.738563820601800e− 03 ∗ z2)+ (4.240508129753556e− 03 ∗ z1) + (1.541011940991888e− 03 ∗ z22)+ (2.115950977442756e− 03 ∗ z1 ∗ z2) + (3.129465795158815e− 04 ∗ z21)+ [−1.316162684352782e− 01, 1.316162684352782e− 01]
The result with symbolic remainder:
(7.879035543566243e− 01) + (−7.425376458676299e− 03 ∗ z2)+ (4.133105242196053e− 03 ∗ z1) + (−7.124243361287374e− 05 ∗ z22)+ (7.930977707887526e− 05 ∗ z1 ∗ z2) + (−2.207266239135132e− 05 ∗ z21)+ [−6.906656051103792e− 02, 6.906656051103792e− 02]

heading: B THEOREM PROOF
text: B.1 PROOF OF SOUNDNESS OF SAMPLING-BASED ERROR ANALYSIS
Proof. The input range of an activation function σj is subdivided into m line segments. Consider the i-th segment [ Y j−Y j m (i− 1) + Y j , Y j−Y j m (i) + Y j ], and let c = Y j−Y j m (i− 1 2 ) + Y j be the center of the segment. The difference between the Bernstein polynomial pjσ and the activation function at the center of the i-th segment is computed as
∣∣pjσ(c)− σj(c)∣∣ . Then, the value of j can be bounded by this difference at the center, as well as the product between the Lipschitz constant of the activation function with respect to this segment Lj and the size of the segment Y j−Y j m , i.e., Lj · Y j−Y j m . The detailed deduction is given below.
|pjσ(x)− σj(x)| =|pjσ(x)− pjσ(c) + pjσ(c)− σj(c) + σj(c)− σj(x)| ≤|pjσ(x)− pjσ(c)|+ |pjσ(c)− σj(c)|+ |σj(c)− σj(x)| Triangle inequality
≤|pjσ(x)− pjσ(c)|+ |pjσ(c)− σj(c)|+ Lj · Y j − Y j
2m Def. of Lipschitz continuity for σj
≤Lj · Y j − Y j
2m + |pjσ(c)− σj(c)|+ Lj · Y j − Y j 2m
Def. of Lipschitz continuity for pjσ [a]
=|pjσ(c)− σj(c)|+ Lj Y j − Y j
m
Note that Bernstein polynomial pjσ has the same Lipschitz constant with σj [a]. Thus we also use Lj to bound |pjσ(x)− pjσ(c)| in the deduction. The error bound over the whole range [Y j , Y j ] should be the largest error bound among all the segments. .
B.2 PROOF OF THEOREM 1
Proof. First, any of our Bernstein overapproximation pσ + Iσ is an input-wise overapproximation for the (component-wise) activation function σ in a layer, i.e., for any input y in the range that Iσ is evaluated, we have that σ(y) ∈ pσ(y) + Iσ. Therefore, the Taylor model (TM) output (pr, Ir) of Algorithm 2 and 3 is a state-wise overapproximation w.r.t. the TM variable z, which is the variable representing the initial state and ranging in the system initial set X0.
We prove Theorem 1 by an induction on the number of control steps j. Assume that N = δc/δ is the number of flowpipes computed in each control step.
Base Case. When j = 0, the TM flowpipes are computed for the reachable set in the first control step and the evolution is under the pure continuous dynamics ẋ = f(x,u0), u̇ = 0 with x(0) ∈ X0 and
u(0) = κ(x(0)). The range of u(0), that is {κ(x(0)) |x(0) ∈ X0}, is overapproximated by a TM (pr(x0), Ir) with z ∈ X0 using Algorithm 2 or 3. Hence, by performing TM flowpipe construction for the ODE ẋ = f(x,u), u̇ = 0 with the initial set x(0) = x0, u(0) ∈ pr(x0) + Ir such that z range in X0, we have that for any i = 1, . . . , N , the i-th TM flowpipe (pi(x0, t), Ii) contains the exact reachable state at the time (i− 1)δ + t. Here, t is the flowpipe local time variable ranging in [0, δ].
Induction. When j > 0, we assume that the local initial set Xj = (p0(x0), I0) is a state-wise overapproximation of the reachable set at the time jδc from any x0 ∈ X0. Then by Algorithm 2 or 3, the obtained TM (pr(x0), Ir) is a state-wise overapproximation for the control input set κ(Xj), i.e., the control output generated by the neural network controller based on the jδc-time state in the execution from any initial state x0 ∈ X0 is contained in the box pr(x0) + Ir. Hence, for any i = 1, . . . , N , the i-th flowpipe (pi(x0, t), Ii) computed for the ODE ẋ = f(x,u), u̇ = 0 with the initial set x(0) ∈ Xj , u(0) ∈ pr(x0) + Ir is a state-wise overapproximation for the reachable set in the time interval of [jδc + (i− 1)δ, jδc + iδ]. More precisely, the reachable state at the time jδc + (i− 1)δ + t for any t ∈ [0, δ] is contained in the box pi(x0, t) + Ii.

heading: C ADDITIONAL EXPERIMENTAL RESULTS
text: Here, we present additional plots of reachable sets computed by different techniques for the benchmarks in Section 4 of the main paper.
For each benchmark, the goal is to check whether the system will reach a given target set. For each tool and in each test, if the computed reachable set overapproximation for the last control step lies entirely in the target set, we consider the tool to have successfully verified the reachability property. If the overapproximation of the reachable set does not intersect with the target set, the tool would have successfully disproved the reachability property. Otherwise, we consider the verification result to be unknown.
The red trajectories are sample system executions and should be contained entirely by the flowpipes computed by each tool. The dark green sets are the flowpipes computed by POLAR. The light green sets are the flowpipes computed by ReachNN* Huang et al. (2019); Fan et al. (2020). The blue sets are the flowpipes computed by Sherlock Dutta et al. (2019). The grey sets are the flowpipes computed by Verisig 2.0 Ivanov et al. (2021b). In some benchmarks, the reachable sets computed by Verisig 2.0 are almost overlapping with the reachable sets computed by POLAR. However, POLAR takes much less time to compute the reachable sets compared to Verisig 2.0 as shown in Table 1 of the main paper. We also show results from NNV Tran et al. (2020) in yellow for some of the benchmarks. For the rest, NNV used up all of the system memory (8GB) and could not finish the computation. Our observations are consistent with those in Ivanov et al. (2021b) where NNV is not able to verify any of these benchmarks. The blue box represents the target set in each test. Except for ex2-sigmoid, POLAR produces the tightest reachable set estimation and successfully proves or disproves the reachability property for all the examples.
D IMPLEMENTATION
The source code can be found in the attached zip. Please follow the instructions described in README.md to run the examples shown in the paper.

Review summaries:
This paper proposes POLAR, a reachability analysis framework for neural network-controlled systems. POLAR improves the previous Taylor model overapproximation approach in two aspects: 1) use Bernstein polynomial to overapproximates activation function; 2) use symbolic remainders to mitigate the so-called wrapping effect, that is, instead of computing an overapproximated interval at each step, keep track of linear transformation matrices to delay overapproximation as late as possible. The experimental evaluation compares POLAR with three other recent work and shows that POLAR achieves the best performance on almost all tasks.

The authors attack a significant problem or bounded-time reachability analysis of neural-network controlled systems (NNCSs). The method allows for verification of decision-making procedures applied by trained neural networks to systems with the known continuous model (continuous plants).  The claimed main novelty of the paper is the application of the basic Taylor Model arithmetic in combination with the univariate Bernstein polynomial interpolation to handle non-differentiable NN activation functions.  The algorithm is tested in a set of benchmarks introduced in the other related work.   The main theoretical result states that the TM flow pipes computed using the POLAR approach are state-wise approximations for the reachable sets.