# PolyARBerNN: A Neural Network Guided Solver and Optimizer for Bounded
Polynomial Inequalities

Wael Fatnassi<EMAIL_ADDRESS>1234-5678-9012 and  Yasser Shoukry
<EMAIL_ADDRESS>0000-0002-8224-8477 University of California,
IrvineCAUSA92697

(2022)

###### Abstract.

Constraints solvers play a significant role in the analysis, synthesis, and
formal verification of complex cyber-physical systems. In this paper, we study
the problem of designing a scalable constraints solver for an important class
of constraints named polynomial constraint inequalities (also known as
nonlinear real arithmetic theory). In this paper, we introduce a solver named
PolyARBerNN that uses convex polynomials as abstractions for highly nonlinears
polynomials. Such abstractions were previously shown to be powerful to prune
the search space and restrict the usage of sound and complete solvers to small
search spaces. Compared with the previous efforts on using convex
abstractions, PolyARBerNN provides three main contributions namely (i) a
neural network guided abstraction refinement procedure that helps selecting
the right abstraction out of a set of pre-defined abstractions, (ii) a
Bernstein polynomial-based search space pruning mechanism that can be used to
compute tight estimates of the polynomial maximum and minimum values which can
be used as an additional abstraction of the polynomials, and (iii) an
optimizer that transforms polynomial objective functions into polynomial
constraints (on the gradient of the objective function) whose solutions are
guaranteed to be close to the global optima. These enhancements together
allowed the PolyARBerNN solver to solve complex instances and scales more
favorably compared to the state-of-art nonlinear real arithmetic solvers while
maintaining the soundness and completeness of the resulting solver. In
particular, our test benches show that PolyARBerNN achieved 100X speedup
compared with Z3 8.9, Yices 2.6, and PVS (a solver that uses Bernstein
expansion to solve multivariate polynomial constraints) on a variety of
standard test benches. Finally, we implemented an optimizer called PolyAROpt
that uses PolyARBerNN to solve constrained polynomial optimization problems.
Numerical results show that PolyAROpt is able to solve high-dimensional and
high order polynomial optimization problems with higher speed compared to the
built-in optimizer in the Z3 8.9 solver.

††copyright: acmcopyright††journalyear: 2022††doi: XXXXXXX.XXXXXXX††journal:
TCPS††journalvolume: 37††journalnumber: 4††article: 1††publicationmonth: 7

##  1\. Introduction

Constraint solvers and optimizers have been used heavily in the design,
synthesis, and verification of cyber-physical systems (Vinod and Oishi, 2018;
mor, 2018; An et al., 2018; Boreale, 2018; Bak et al., 2019; Xue et al., 2019;
Dutta et al., 2019; Marcucci and Tedrake, 2019; Kochdumper and Althoff, 2020;
Sadraddini and Tedrake, 2020; Rabi, 2020). Examples includes verification of
neural network controlled autonomous systems (Sun et al., 2019), formal
verification of human-robot interaction in healthcare scenarios (Lestingi et
al., 2020), automated synthesis for distributed cyber-physical systems (Roy et
al., 2018), design for cyber-physical systems under sensor attacks (Shoukry et
al., 2018), air traffic management of unmanned aircraft systems (Munoz, 2015),
software verification for the next generation space-shuttle (Nelson and
Pecheur, 2002), and conflict detection for aircraft (Narkawicz and Munoz,
2012).

In this paper, we will focus on the class of general multivariate polynomial
constraints (also known as nonlinear real arithmetic). Multivariate polynomial
constraints appear naturally in the design, synthesis, and verification of
these systems. It is not then surprising that the amount of attention given to
this problem in the last decade, as evidenced by the amount of off-the-self
solvers that are designed to solve feasibility and optimization problems over
general multivariate polynomial constraints, including Z3 (De Moura and
Bjørner, 2008), Coq (Mahboubi, 2006), Yices (Dutertre, 2014), PVS (Munoz and
Narkawicz, 2013), Cplex, (Manual, 1987), CVXOPT (Vandenberghe, 2010), and
Quadprog (Goldfarb and Idnani, 1983). Regardless of their prevalence in
several synthesis and verification problems, well-known algorithms—that are
capable of solving a set of polynomial constraints—are shown to be doubly
exponential (England and Davenport, 2016), placing a significant challenge to
design efficient solvers for such problems.

Recently, neural networks (NNs) have shown impressive empirical results in
approximating unknown functions. This observation motivated several
researchers to ask how to use NNs to tame the complexity of NP-hard problems.
Examples are the use of NNs to design scalable solvers for program synthesis
(Devlin et al., 2017), traveling salesman problem (Bello et al., 2016), first-
order theorem proving (Irving et al., 2016), higher-order theorem proving
(Kaliszyk et al., 2017), and Boolean satisfiability (SAT) problems (Selsam et
al., 2019). While several of these solvers sacrifice either soundness or
correctness guarantees, we are interested in this paper on using such
empirically powerful NNs to design a sound and complete solver for nonlinear
real arithmetic.

In addition to NNs, polynomials constitute a rich class of functions for which
several approximators have been studied. Two of the most famous approximators
for polynomials are Taylor approximation and Bernstein polynomials. These two
approximators have been successfully used in solvers like Coq and PVS
(Brisebarre et al., 2012; Munoz and Narkawicz, 2013). This opens the question
of how to combine all those approximation techniques, i.e., NNs, Taylor, and
Bernstein approximations, to come up with a scalable solver that can reason
about general multivariate polynomial constraints. We introduce PolyARBerNN, a
novel sound and complete solver for polynomial constraints that combines these
three function approximators (NNs, Taylor, and Bernstein) to prune the search
space and produce small enough instances in which existing sound and complete
solvers (based on the well-known Cylindrical Algebraic Decomposition
algorithm) can easily reason about. In general, we provide the following
contributions:

  * •

We introduced a novel NN-guided abstraction refinement process in which a NN
is used to guide the use of Taylor approximations to find a solution or prune
the search space. We analyzed the theoretical characteristics of such a NN and
provided empirical evidence on the generalizability of the trained NN in terms
of its ability to guide the abstraction refinement process for unseen
polynomials with various numbers of variables and orders.

  * •

We complement the NN-guided abstraction refinement with a state-space pruning
phase using Bernstein approximations that accelerates the process of removing
portions of the state space in which the sign of the polynomial does not
change.

  * •

We validated our approach by first comparing the scalability of the proposed
PolyARBerNN solver with respect to PVS, a library that uses Bernstein
expansion to solve polynomial constraints. Second, we compared the execution
times of the proposed tool with the latest versions of the state-of-the-art
nonlinear arithmetic solvers, such as Z3 8.9, Yices 2.6 by varying the order,
the number of variables, and the number of the polynomial constraints for
instances when a solution exists and when a solution does not exist. We also
compared the scalability of the solver against Z3 8.9 and Yices 2.9 on the
problem of synthesizing a controller for a cyber-physical system.

  * •

We proposed PolyAROpt, an optimizer that uses PolyARBerNN to solve constrained
multivariate polynomial optimization problems. Our theoretical analysis shows
that PolyAROpt is capable of providing solutions that are $\epsilon$ close to
the global optima (for any $\epsilon>0$ chosen by the user). Numerical results
show that PolyAROpt solves high-dimensional and high-order optimization
problems with high speed compared to the built-in optimizer in Z3 8.9 solver.
We also validated the effectiveness of PolyAROpt on the problem of computing
the reachable sets of polynomial dynamical systems.

Related work: Cylindrical algebraic decomposition (CAD) was introduced by
Collins (Collins, 1975) in 1975 and is considered to be the first algorithm to
effectively solve general polynomial inequality constraints. Several
improvements were introduced across the years to reduce the high time
complexity of the CAD algorithm (Hong, 1990; McCallum, 1998; Brown, 2001).
Although the CAD algorithm is sound and complete, it scales poorly with the
number of polynomial constraints and their order. Other techniques to solve
general polynomial inequality constraints include the use of transformations
and approximations to scale the computations. For instance, the authors in
(Munoz and Narkawicz, 2013) incorporated Bernstein polynomials in the
Prototype Verification System (PVS) theorem prover; these developments are
publicly available in the NASA PVS Library. The library uses the range
enclosure propriety of Bernstein polynomials to solve quantified polynomial
inequalities. However, the library is not complete for non-strict inequalities
(Munoz and Narkawicz, 2013) and is not practical for higher dimensional
polynomials. Another line of work that is related to our work is the use of
machine learning to solve combinatorial problems (Gasse et al., 2019; Selsam
et al., 2019). In particular, the authors in (Gasse et al., 2019) proposed a
graph convolutional neural network (GCNN) to learn heuristics that can
accelerate mixed-integer linear programming (MILP) solvers. Similarly, the
NeuroSAT solver (Selsam et al., 2019) uses a message-passing neural network
(MPNN) to solve Boolean SAT problems. The authors of (Selsam et al., 2019)
showed that NeuroSAT generalizes to novel distributions after training only on
random SAT problems. Nevertheless, NeuroSAT is not competitive with state-of-
art SAT solvers and it does not have a correctness guarantee.

##  2\. Problem Formulation

Notation: We use the symbols $\mathbb{N}$ and $\mathbb{R}$ to denote the set
of natural and real numbers, respectively. We denote by
$x=\big{(}x_{1},x_{2},\cdots,x_{n}\big{)}\in\mathbb{R}^{n}$ the vector of
real-valued variables, where $x_{i}\in\mathbb{R}$. We denote by
$I_{n}(\underline{d},\overline{d})=\big{[}\underline{d}_{1},\overline{d}_{1}\big{]}\times\cdots\times$
$\big{[}\underline{d}_{n},\overline{d}_{n}\big{]}\subset\mathbb{R}^{n}$ the
$n$-dimensional hyperrectangle where
$\underline{d}=\left(\underline{d}_{1},\cdots,\underline{d}_{n}\right)$ and
$\overline{d}=\left(\overline{d}_{1},\cdots,\overline{d}_{n}\right)$ are the
lower and upper bounds of the hyperrectangle, respectively. For a real-valued
vector $x=\big{(}x_{1},x_{2},\cdots,x_{n}\big{)}\in\mathbb{R}^{n}$ and an
index-vector $K=\left(k_{1},\cdots,k_{n}\right)\in\mathbb{N}^{n}$, we denote
by $x^{K}\in\mathbb{R}$ the scalar $x^{K}=x_{1}^{k_{1}}\cdots x_{n}^{k_{n}}$.
Given two multi-indices $K=\left(k_{1},\cdots,k_{n}\right)\in\mathbb{N}^{n}$
and $L=\left(l_{1},\cdots,l_{n}\right)\in\mathbb{N}^{n}$, we use the following
notation throughout this paper:
$K+L=\left(k_{1}+l_{1},\cdots,k_{n}+l_{n}\right)$, ${L\choose K}={l_{1}\choose
k_{1}}\cdots{l_{n}\choose k_{n}}$, and $\sum\limits_{K\leq
L}=\sum\limits_{k_{1}\leq l_{1}}\cdots\sum\limits_{k_{n}\leq l_{n}}$. A real-
valued multivariate polynomial $p:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is
defined as:

$\displaystyle p(x_{1},\ldots,x_{n})$
$\displaystyle\;=\;\sum_{k_{1}=0}^{l_{1}}\sum_{k_{2}=0}^{l_{2}}\ldots\sum_{k_{n}=0}^{l_{n}}a_{(k_{1},\ldots,k_{n})}x_{1}^{k_{1}}x_{2}^{k_{2}}\ldots
x_{n}^{k_{n}}\;=\;\sum\limits_{K\leq L}a_{K}x^{K},$

where $L=(l_{1},l_{2},\ldots,l_{n})$ is the maximum degree of $x_{i}$ for all
$i=1,\ldots,n$. We denote by
$a_{p}=(a_{(0,0,\ldots,0)},\ldots,a_{(l_{1},l_{2},\ldots,l_{n})})$ the vector
of all the coefficients of polynomial $p$. We denote the space of multivariate
polynomials with coefficients in $\mathbb{R}$ by
$\mathbb{R}[x_{1},x_{2},\cdots,x_{n}]$. Given a real-valued function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, we denote by
$L^{-}_{0}\left(f\right)$ and $L^{+}_{0}\left(f\right)$ the zero sublevel and
zero superlevel sets of f, i.e.,:

$\displaystyle
L^{-}_{0}\left(f\right)=\\{x\in\mathbb{R}^{n}\big{|}f\big{(}x\big{)}\leq
0\\},\qquad
L^{+}_{0}\left(f\right)=\\{x\in\mathbb{R}^{n}\big{|}f\big{(}x\big{)}\geq
0\\}.$

Finally, a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is called
Lipschitz continuous if there exists a positive real constant $\omega_{f}\geq
0$ such that, for all $x_{1}\in\mathbb{R}^{n}$ and $x_{2}\in\mathbb{R}^{n}$,
the following holds:

$\left\lVert f(x_{1})-f(x_{2})\right\rVert\leq\omega_{f}\left\lVert
x_{1}-x_{2}\right\rVert$

Main Problem: In this paper, we focus on two problems namely (Problem 1) the
feasibility problems that involve multiple polynomial inequality constraints
with input ranges confined within closed hyperrectangles and (Problem 2) the
constrained optimization problem which aims to maximize (or minimize) a
polynomial objective function subject to other polynomial inequality
constraints and input range constraints.

###### Problem 1.

$\displaystyle\exists x\in I_{n}(\underline{d},\overline{d})\quad\text{such
that:}\quad$ $\displaystyle
p_{1}\left(x_{1},\cdots,x_{n}\right)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 0\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
$\displaystyle\qquad\qquad\vdots$ $\displaystyle
p_{m}\left(x_{1},\cdots,x_{n}\right)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 0$

where
$p_{i}\left(x\right)=p_{i}\left(x_{1},\cdots,x_{n}\right)\in\mathbb{R}[x_{1},x_{2},\cdots,x_{n}]$
is a polynomial over variables $x_{1},\cdots,x_{n}$. Without loss of
generality, $p_{i}\big{(}x\big{)}$ $\geq\leavevmode\nobreak\ 0$ and
$p_{i}\left(x\right)\leavevmode\nobreak\ =\leavevmode\nobreak\ 0$ can be
encoded using the constraints above. Similarly, given a polynomial objective
function $p(x)\in\mathbb{R}[x_{1},x_{2},\cdots,x_{n}]$, we define the
optimization problem as:

###### Problem 2.

$\displaystyle\min\limits_{x\in
I_{n}(\underline{d},\overline{d})}\leavevmode\nobreak\ $ $\displaystyle
p\left(x\right)\leavevmode\nobreak\ \leavevmode\nobreak\
[\textbf{or}\leavevmode\nobreak\ \max\limits_{x\in
I_{n}(\underline{d},\overline{d})}\leavevmode\nobreak\
p\left(x\right)]\qquad\qquad\qquad$ subject to: $\displaystyle
p_{1}\left(x_{1},\cdots,x_{n}\right)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 0,$ $\displaystyle\qquad\qquad\vdots$ $\displaystyle
p_{m}\left(x_{1},\cdots,x_{n}\right)\leavevmode\nobreak\
\leq\leavevmode\nobreak\ 0$

##  3\. Convex Abstraction Refinement: Benefits and Drawbacks

In this section, we overview our previously reported framework for using
convex abstraction refinement process introduced in (Yasser Shoukry, 2021)
along with some drawbacks that motivate the need for the proposed framework.

###  3.1. Overview of Convex Abstraction Refinement

Sound and complete algorithms that solve Problem 1 are known to be doubly
exponential in $n$ with a total running time that it is bounded by
$\left(m\leavevmode\nobreak\ \overline{deg}\right)^{2^{n}}$ (England and
Davenport, 2016), where $\overline{deg}$ is the maximum degree among the
polynomials $p_{1},\ldots,p_{m}$. Since the complexity of the problem grows
exponentially, it is useful to remove (or prune) subsets of the search space
in which the solution is guaranteed not to exist. Since Problem 1 asks for an
$x$ in $\mathbb{R}^{n}$ for which all the polynomials are negative, a solution
does not exist in subsets of $\mathbb{R}^{n}$ at which one of the polynomials
$p_{i}$ is always positive (i.e., $L^{+}_{0}\left(p_{i}\right)$). In the same
way, finding regions of the input space for which some of the polynomials are
negative $L^{-}_{0}\left(p_{i}\right)$ helps with finding the solution faster.

To find subsets of $L^{+}_{0}\left(p_{i}\right)$ and
$L^{-}_{0}\left(p_{i}\right)$ efficiently, the use of “convex abstractions” of
the polynomials was previously proposed by the authors in (Yasser Shoukry,
2021). Starting from a polynomial $p_{i}(x)\in\mathbb{R}[x]$ and a
hyperrectangle $I_{n}\subset\mathbb{R}^{n}$, the framework in (Yasser Shoukry,
2021) computes two quadratic polynomials $O^{p_{i}}_{j}$ and $U^{p_{i}}_{j}$
such that:

(1) $\displaystyle U^{p_{i}}_{j}(x)\leq p(x)\leq
O^{p_{i}}_{j}(x),\qquad\forall x\in I_{n},$

where $O$ and $U$ stands for Over-approximate and Under-approximate quadratic
polynomials, respectively, and the subscript $j$ in $O^{p_{i}}_{j}(x)$ and
$U^{p_{i}}_{j}(x)$ encodes the iteration index of the abstraction refinement
process. It is easy to notice that the zero superlevel set of
$U^{p_{i}}_{j}(x)$ is a subset of $L^{+}_{0}(p_{i})$, i.e.,
$L^{+}_{0}(U^{p_{i}}_{j})\subseteq L^{+}_{0}(p_{i})$. Similarly, the zero
sublevel set of $O^{p_{i}}_{j}(x)$ is a subset of $L^{-}_{0}(p_{i})$, i.e.,
$L^{-}_{0}(O^{p_{i}}_{j})\subseteq L^{-}_{0}(p_{i})$. Moreover, being convex
polynomials, identifying the zero superlevel sets and zero sublevel sets of
$O^{p_{i}}_{j}(x)$ and $U^{p_{i}}_{j}(x)$ can be computed efficiently using
convex programming tools. By iteratively refining these upper and lower convex
approximations, the framework in (Yasser Shoukry, 2021) was able to rapidly
prune the search space until regions with relatively small volumes are
identified, at which sound and complete tools such as Z3 8.9 and Yices 2.6
(which are based on the Cylindrical Algebraic Decomposition algorithm) are
used to search these small regions, efficiently, to find a solution. It is
important to notice that these solvers (especially Yices) are optimized for
the cases when the search space is a bounded hyperrectangle.

Figure 1. Exemplary cases where abstracting higher order polynomial (black
curves) using convex approximations fails to provide helpful information: Top-
Left: under-approximation (green curve) is entirely negative and hence fails
to identify any subsets of $L^{+}_{0}\left(p\right)$. Top-Right: over-
approximation (red curve) is entirely positive and hence fails to identify
subsets of $L^{-}_{0}\left(p\right)$. Bottom: under/over approximations failed
to identify polynomials that are consistently positive (left) or negative
(right).

###  3.2. Drawbacks of Convex Abstraction Refinement

Although the prescribed convex abstraction refinement process was shown to
provide several orders of magnitude speedup compared to the state-of-the-art
(Yasser Shoukry, 2021), it adds unnecessary overhead in certain situations. In
particular, and as shown in Figure 1, the quadratic abstractions
$O^{p_{i}}_{j}(x)$ and $U^{p_{i}}_{j}(x)$ may fail to identify meaningful
subsets of $L^{-}_{0}\left(p_{i}\right)$ and $L^{+}_{0}\left(p_{i}\right)$.
One needs to split the input region to tighten the over-/under-approximation
in such cases. Indeed, applying the convex abstraction refinement process,
above, may lead to several unnecessary over-approximations or under-
approximations until a tight one that prunes the search space is found. These
drawbacks call for a methodology that is capable of:

  1. (1)

Guiding the abstraction refinement process: To reduce the number of
unnecessary computations of over/under approximations, one needs a heuristic
that guides the convex abstraction refinement process. In particular, such a
heuristic needs to consider the properties of the polynomials and the input
region to estimate the volume of the sets that the convex under/over-
approximation will identify.

  2. (2)

Alternative Abstraction: As shown in Figure 1 (bottom), abstracting high-order
polynomials using convex ones may fail to identify easy cases when the
polynomial is strictly positive or negative. Therefore, it is beneficial to
use alternative ways to abstract high-order polynomials that can augment the
convex abstractions.

Designing a strategy that addresses the two requirements above is the main
topic for the following two sections.

##  4\. Neural Network Guided Convex Abstraction Refinement

In this section, we are interested in designing and training a Neural Network
(NN) that can be used to guide the abstraction refinement process. Such NN can
be used as an oracle by the solver to estimate the volume of the zero
super/sub-level sets (for each polynomial) within a given region
$I_{n}(\underline{d},\overline{d})$ and select the best approximation strategy
out of three possibilities namely: (i) apply convex under-approximation, (ii)
apply convex over-approximation, and (iii) split the region to allow for finer
approximations in the subsequent iterations of the solver. In this section, we
aim to develop a scientific methodology that can guide the design of such NN.

###  4.1. On the relation between the NN architecture and the characteristics
of the polynomials:

In this subsection, we aim to understand how the properties of the polynomials
affect the design of the NN. We start by reviewing the following result from
the machine learning literature:

###### Theorem 1 (Theorem 1.1 (Shen et al., 2019)).

There exists a Rectifier Linear Unit (ReLU)-based neural network $\phi$ that
can estimate a continuous function $f$ such that the estimation error is
bounded by:

$||\phi-f||\leq\omega_{f}\;\sqrt{d}\;O(N^{-2/d}L^{-2/d})$

where $N,L,d$ are the neural network depth, the neural network width, and the
number of neural network inputs, respectively, and $\omega_{f}$ is the
Lipschitz constant of the function $f$. Moreover, this bound is nearly tight.

The above result can be interpreted as follows. The depth $N$ and width $L$ of
a neural network depend on the rate of change of the underlying function
(captured by its Lipschitz constant $\omega_{f}$). That is, if we use a NN to
estimate a function with a high $\omega_{f}$, then one needs to increase the
depth $N$ and width $L$ of the NN to achieve an acceptable estimation error.

Now we aim to connect the result above with the characteristics of the
polynomials. To that end, we recall the definition of “condition numbers” of a
polynomial (Farouki and Rajan, 1987):

###### Definition 0.

Given a polynomial $p\left(x\right)=\sum\limits_{K}a_{K}x^{K}$ and a root
$x_{0}$ of $p$, the quantity $C_{a_{p}}\left(x_{0}\right)$ is called the
condition number for the root $x_{0}$. The condition number characterizes the
sensitivity of the root $x_{0}$ to a perturbation of the coefficients $a_{p}$.
That is, if we allow a random perturbation of a fixed relative magnitude
$\epsilon=\lvert\frac{\delta a_{K}}{a_{K}}\rvert$ in each coefficient $a_{K}$
in $a_{p}$, then the magnitude of the maximum displacement $\delta x_{0}$ of a
root $x_{0}$ is bounded as: $\lvert\delta x_{0}\rvert\leq
C_{a_{p}}\left(x_{0}\right)\epsilon$. For a polynomial with multiple roots,
then we define the condition number of the polynomial $\overline{C}_{a_{p}}$
as the largest $C_{a_{p}}\left(x_{0}\right)$ among all roots, i.e.,
$\overline{C}_{a_{p}}=sup_{x_{0}\in\\{x|p(x)=0\\}}C_{a_{p}}\left(x_{0}\right)$.

We are now ready to present our first theoretical result that connects the
condition number of polynomials to the NN architecture. As stated before, we
are interested in designing an NN that can estimate the zero sub/super level
volume set within a given region. We show that the larger the condition
number, the larger the neural network depth and width, as captured by the
following result.

###### Theorem 3.

Given a polynomial $p$ with coefficients $a_{p}$ and a region
$I_{n}(\underline{d},\overline{d})$. There exists a neural network
$NN(a_{p},I_{n})$ that estimates the volume of zero sub/super level sets from
the polynomial coefficients $a_{p}$. The Lipschitz constant of this
$NN(a_{p},I_{n})$, denoted by $\omega_{NN}$ is bounded by
$\mathcal{O}(n_{r}\overline{n}_{r}\overline{C}_{a_{p}})$ where
$\overline{C}_{a_{p}}$ is the condition number of the polynomial $p$,
$\overline{n}_{r}=max(n,n_{r})$ and $n_{r}$ is the number of roots of the
polynomial $p$.

To prove the result, we will proceed with an existential argument. We will
show that a NN that matches the properties above exists without constructing
such a NN. As shown in Figure 2, the neural network $NN(a_{p},I_{n})$ consists
of multiple sub-neural networks. In particular, the first sub-neural network
$NN_{a_{p}\rightarrow X_{0}}$ computes all the roots
$X_{0}=(x_{0}^{1},\ldots,x_{0}^{n_{r}})$ of the polynomial (where $n_{r}$ is
the number of roots) from the coefficients $a_{p}$, i.e.:

(2) $\displaystyle X_{0}=NN_{a_{p}\rightarrow X_{0}}(a_{p}).$

Note that $NN_{a_{p}\rightarrow X_{0}}$ does not depend on the region $I_{n}$
and hence the roots $X_{0}$ may not lie inside the region $I_{n}$. Moreover,
Theorem 3 asks for a NN that estimates the volume of the zero sub/super level
sets and not the location of the roots. To that end, our strategy is to split
the region $I_{n}$ into sub-regions of fixed volume and check if a root lies
within each of these sub-regions. If a sub-region does not have a root (i.e.,
there is no zero crossing inside this sub-region) and the evaluation of the
polynomial at any point in this region turns to be positive, then this sub-
region belongs to the super level set of $p$ and similarly for the sublevel
set of $p$. By counting the number of the sub-regions with no zero crossings
and multiplying this count by the volume of these sub-regions, we can provide
an estimate of the sub/super level sets. Such a process can be performed using
the following three sub-neural networks:

  * •

The sub-neural network $NN_{I_{n}\rightarrow I^{i}_{n}}$ splits the region
$I_{n}$ into $l$ sub-regions $I^{1}_{n},\ldots,I^{l}_{n}$ and return the
bounds of the $i$th sub-region, i.e.:

(3) $\displaystyle(\underline{d}^{i},\overline{d}^{i})$
$\displaystyle=NN_{I_{n}\rightarrow I^{i}_{n}}(I_{n}),\qquad
i\in\\{1,\ldots,l\\}.$

  * •

The sub-neural network $NN_{X_{0}\rightarrow ZC_{i}}$ checks the location of
the roots $(x^{1}_{0},\ldots,x_{0}^{n_{r}})$ and returns a binary indicator
variable $ZC_{i}$ that indicates whether a zero-crossing takes place within
the $i$th sub-region or whether the polynomial is always positive/negative
within the $i$th sub-region, i.e.:

(4) $\displaystyle ZC_{i}(a_{p},I^{i}_{n})$
$\displaystyle=NN_{X_{0}\rightarrow ZC_{i}}\left(NN_{a_{p}\rightarrow
X_{0}}(a_{p}),NN_{I_{n}\rightarrow I^{i}_{n}}(I_{n})\right).$

  * •

The final output $NN(a_{p},I_{n})$ is computed using the sub-neural network
$NN_{ZC\rightarrow L^{+}/L^{-}}$ which counts the number of regions that has
no zero-crossing (using the indicators $ZC_{1},\ldots,ZC_{l}$) and compute the
estimate of the zero sub/super level sets, i.e.:

(5) $\displaystyle NN(a_{p},I_{n})$ $\displaystyle=NN_{ZC\rightarrow
L^{+}/L^{-}}\left(ZC_{1}(a_{p},I^{1}_{n}),\ldots,ZC_{l}(a_{p},I^{l}_{n})\right).$

The Lipschitz constants of these sub-neural networks are captured by the
following four propositions whose proof can be found in the appendix.

###### Proposition 0.

Consider the sub-neural network $NN_{a_{p}\rightarrow X_{0}}(a_{p})$ defined
in (2). The Lipschitz constant of $NN_{a_{p}\rightarrow X_{0}}(a_{p})$ is
bounded by $\mathcal{O}(n_{r}\overline{C}_{a_{p}})$.

###### Proposition 0.

Consider the sub-neural network $NN_{I_{n}\rightarrow I^{i}_{n}}$ defined in
(3). The Lipschitz constant of $NN_{I_{n}\rightarrow I^{i}_{n}}$ is bounded by
$\mathcal{O}(n)$.

Figure 2. The architecture of the neural network $NN(a_{p},I_{n})$ used to
prove the Theorem 3.

###### Proposition 0.

Consider the sub-neural network $NN_{X_{0}\rightarrow
ZC_{i}}\left(X_{0},I^{i}_{n}\right)$ defined in (4). The Lipschitz constant of
$NN_{X_{0}\rightarrow ZC_{i}}\left(X_{0},I^{i}_{n}\right)$ is bounded by
$\mathcal{O}(n_{r})$.

###### Proposition 0.

Consider the sub-neural network $NN_{ZC\rightarrow L^{+}/L^{-}}$ defined in
(5). The Lipschitz constant of $NN_{ZC\rightarrow L^{+}/L^{-}}$ is bounded by
$\mathcal{O}(1)$.

###### Proof of Theorem 3.

Consider the NN shown in Figure 2 and defined using equations (2)-(5). To
bound the Lipschitz constant of $NN(a_{p},I_{n})$, we consider two sets of
inputs $(a_{p},I_{n})$ and $(a^{\prime}_{p},I^{\prime}_{n})$ as follows:

$\displaystyle\left\lVert
NN(a^{\prime}_{p},I^{\prime}_{n})-NN(a_{p},I_{n})\right\rVert_{2}$
$\displaystyle=\left\|NN_{ZC\rightarrow
L^{+}/L^{-}}\left(ZC_{1}(a^{\prime}_{p},I^{\prime
1}_{n}),\ldots,ZC_{l}(a^{\prime}_{p},I^{\prime l}_{n})\right)\right.$ (6)
$\displaystyle\left.\qquad\qquad-NN_{ZC\rightarrow
L^{+}/L^{-}}\left(ZC_{1}(a_{p},I^{1}_{n}),\ldots,ZC_{l}(a_{p},I^{l}_{n})\right)\right\|_{2},$
(7)
$\displaystyle\leq\mathcal{O}(1)\left\lVert\left(ZC_{1}(a^{\prime}_{p},I^{\prime
1}_{n})-ZC_{1}(a_{p},I^{1}_{n}),\cdots,ZC_{l}(a^{\prime}_{p},I^{\prime
l}_{n})-ZC_{l}(a_{p},I^{l}_{n})\right)\right\rVert_{2},$ (8)
$\displaystyle=\mathcal{O}(1)\bigg{(}\sum\limits_{i=1}^{l}\left\lVert
ZC_{i}(a^{\prime}_{p},I^{\prime
i}_{n})-ZC_{i}(a_{p},I^{i}_{n})\right\rVert^{2}_{2}\bigg{)}^{\frac{1}{2}}.$

where (7) follows from Proposition 7. Now, we upper bound $\left\lVert
ZC_{i}(a^{\prime}_{p},I^{\prime
i}_{n})-ZC_{i}(a_{p},I^{i}_{n})\right\rVert_{2}$ as follows:

$\displaystyle\left\lVert ZC_{i}(a^{\prime}_{p},I^{\prime
i}_{n})-ZC_{i}(a_{p},I^{i}_{n})\right\rVert^{2}_{2}$
$\displaystyle=\left\|NN_{X_{0}\rightarrow ZC_{i}}\left(NN_{a_{p}\rightarrow
X_{0}}(a^{\prime}_{p}),NN_{I_{n}\rightarrow
I^{i}_{n}}(I^{\prime}_{n})\right)\right.$ (9)
$\displaystyle\left.\qquad\qquad\qquad-NN_{X_{0}\rightarrow
ZC_{i}}\left(NN_{a_{p}\rightarrow X_{0}}(a_{p}),NN_{I_{n}\rightarrow
I^{i}_{n}}(I_{n})\right)\right\|^{2}_{2}$ (10)
$\displaystyle\leq\mathcal{O}(n_{r})\left\lVert\left(NN_{a_{p}\rightarrow
X_{0}}(a^{\prime}_{p})-NN_{a_{p}\rightarrow X_{0}}(a_{p}),NN_{I_{n}\rightarrow
I^{i}_{n}}(I^{\prime}_{n})-NN_{I_{n}\rightarrow
I^{i}_{n}}(I_{n})\right)\right\rVert^{2}_{2},$
$\displaystyle=\mathcal{O}(n_{r})\left(\left\lVert\left(NN_{a_{p}\rightarrow
X_{0}}(a^{\prime}_{p})-NN_{a_{p}\rightarrow
X_{0}}(a_{p})\right)\right\rVert^{2}_{2}+\left\lVert NN_{I_{n}\rightarrow
I^{i}_{n}}(I^{\prime}_{n})-NN_{I_{n}\rightarrow
I^{i}_{n}}(I_{n})\right\rVert^{2}_{2}\right),$ (11)
$\displaystyle\leq\mathcal{O}(n_{r})\left(\mathcal{O}(n_{r}\overline{C}_{a_{p}})\left\lVert
a^{\prime}_{p}-a_{p}\right\rVert^{2}_{2}+\mathcal{O}(n)\left\lVert
I^{\prime}_{n}-I_{n}\right\rVert^{2}_{2}\right)$ (12)
$\displaystyle=\mathcal{O}(n_{r}\overline{n}_{r}\overline{C}_{a_{p}})\left\lVert(a^{\prime}_{p},I^{\prime}_{n})-(a_{p},I_{n})\right\rVert_{2}^{2}.$

where (10) follows from Proposition 6; (11) follows from Propositions 5 and 4
along with the definition of $\overline{n}_{r}=\max(n,n_{r})$. Substituting
(12) in (8) and noticing that $l$ is a constant that does not depend on $n$
yields:

(13) $\displaystyle\left\lVert
NN(a^{\prime}_{p},I^{\prime}_{n})-NN(a_{p},I_{n})\right\rVert_{2}\leq\mathcal{O}(n_{r}\overline{n}_{r}\overline{C}_{a_{p}})\left\lVert(a^{\prime}_{p},I^{\prime}_{n})-(a_{p},I_{n})\right\rVert_{2},$

from which we conclude that the Lipschitz constant of $NN(a_{p},I_{n})$ is in
the order of $\mathcal{O}(n_{r}\overline{n}_{r}\overline{C}_{a_{p}})$ which
concludes the proof of Theorem 3.

∎

It follows from Theorem 1 and Theorem 3 that the higher the condition number
of a polynomial $\overline{C}_{a_{p}}$, the larger the network width and depth
needed to estimate the volume of zero sub/super level sets with high accuracy.
Unfortunately, the power basis representation of polynomials (i.e.,
representing the polynomial as a summation $\sum_{K\leq L}a_{K}x^{K}$) is
shown to be an unstable representation with extremely large condition numbers
(Farouki and Rajan, 1987) which may necessitate neural networks with
substantial architecture.

###  4.2. Bernstein Polynomials: A Robust Representation of Polynomials

Motivated by the challenge above, we seek a representation of polynomials that
is more robust to changes in coefficients, i.e., we seek a representation in
which the roots of the polynomial change slowly with changes in the
coefficients (and hence smaller condition numbers $\overline{C}_{a_{p}}$ and a
smaller NN to estimate the volume of the sub/super level sets). We start with
the following definition.

###### Definition 0.

Let $p\left(x\right)=\sum\limits_{K\leq
L}a_{K}x^{K}\in\mathbb{R}[x_{1},\ldots,x_{n}]$ be a multivariate polynomial
over a hyperrectangle $I_{n}(\underline{d},\overline{d})$ and of a maximal
degree $L=\left(l_{1},\cdots,l_{n}\right)\in\mathbb{N}^{n}$. The polynomial:

(14) $\displaystyle B_{p,L}\left(x\right)$ $\displaystyle=\sum\limits_{K\leq
L}b_{K,L}Ber_{K,L}\left(x\right),$

is called the Bernstein polynomial of $p$, where $Ber_{K,L}\left(x\right)$ and
$b_{K,L}$ are called the Bernstein basis and Bernstein coefficients of $p$,
respectively, and are defined as follows:

(15) $\displaystyle Ber_{K,L}\left(x\right)$ $\displaystyle={L\choose
K}x^{K}\left(1-x\right)^{L-K},\qquad
b_{K,L}=\sum\limits_{J=(0,\ldots,0)}^{K}\frac{{K\choose J}}{{L\choose
J}}\left(\overline{d}-\underline{d}\right)^{J}\sum\limits_{I=J}^{L}{I\choose
J}\underline{d}^{I-J}a_{I}.$

The Bernstein representation is known to be the most robust representation of
polynomials which is captured by the next result (Farouki and Rajan, 1987).

###### Theorem 9 (Theorem (Farouki and Rajan, 1987)).

The Bernstein basis is optimally stable, i.e. there exists no other basis with
a condition number smaller than the condition number of the Bernstein
coefficients $\overline{C}_{b_{p}}$, where
$b_{p}=(b_{(0,0,\ldots,0),L},\ldots,b_{(l_{1},l_{2},\ldots,l_{n}),L})$ is the
vector of all the Bernstein coefficients of polynomial $p$.

Theorems 1-9 point to the optimal way of designing the targeted neural
network. Such a neural network needs to take as input the Bernstein
coefficients $b_{p}$ instead of the power basis coefficients $a_{p}$. To
validate this conclusion, we report empirical evidence in Table 1. In this
numerical experiment, we trained two neural networks with the same exact
architecture, using the same exact number of data points, and both networks
have the same number of inputs. Both neural networks are trained to estimate
whether a zero-crossing occurs in a region (recall from our analysis in
Theorem 3 that the Lipschitz constant of this NN is equal to the condition
number of the polynomial). The only difference is that one neural network is
trained using power basis coefficients $a_{p}$ (column 3 of Table 1) while the
second is trained using Bernstein basis coefficients $b_{p}$ (column 5 of
Table 1). The coefficients are randomly generated via a uniform distribution
between $-0.1$ and $0.1$, i.e., $\mathcal{U}\left(-0.1,0.1\right)$. We
generated 40000 training samples and 10000 validation samples for both bases.
We evaluate the trained NN on three different benchmarks for the two bases.
Each evaluation benchmark has 10000 samples. The results are summarized in
Table 1. As it can be seen from Table 1, the NN trained with Bernstein
coefficients generalizes better than the NN trained with power basis
coefficients as reflected by the empirical “Accuracy” during evaluation. This
empirical evidence matches our analysis in Theorem 3 along with the insights
of Theorem 1 and Theorem 9.

Table 1. Evaluation of three trained neural networks on three different
benchmarks for the different polynomial basis. Each benchmark has 10000
samples. The coefficients of the polynomial within each basis are generated
following a uniform distribution given in the table.

Benchmark Coefficients Power Basis Bernstein Basis Reduced Bernstein Basis
Accuracy Overhead Accuracy Overhead Accuracy Overhead 1
$\mathcal{U}\left(-0.1,0.1\right)$ $46\%$ 0 [s] $91\%$ 0.01 [s] $82\%$ 0.002
[s] 2 $\mathcal{U}\left(-0.5,0.5\right)$ $32\%$ 0 [s] $87\%$ 0.03 [s] $79\%$
0.005 [s] 3 $\mathcal{U}\left(-1,1\right)$ $30\%$ 0 [s] $88\%$ 0.04 [s] $80\%$
0.007 [s]

##  5\. Taming the Complexity of Computing Bernstein Coefficients

In section 4, we concluded that Bernstein’s representation has a smaller
condition compared to other representations, which helps build a more
efficient NN. Nevertheless, computing this representation adds a significant
overhead even by using the most efficient algorithms to calculate these
coefficients (Smith, 2009; Ray and Nataraj, 2012). For example, computing all
the Bernstein coefficients of a $6^{th}$-dimensional polynomial with $7^{th}$
order using Matrix method and Garloff’s methods (Smith, 2009; Ray and Nataraj,
2012) require $1.1e07$ and $7.1e06$ summation and multiplication operations
(Ray and Nataraj, 2012). To exacerbate the problem, the Bernstein coefficients
depend on the region $I_{n}$ and need to be recomputed in every iteration of
the abstraction refinement process. Reducing such overhead is the main focus
of this section.

###  5.1. Range Enclosure Property of Bernstein polynomials

Given a multivariate polynomial $p\left(x\right)$ that is defined over the
$n$-dimensional box $I_{n}(\underline{d},\overline{d})$, we can bound the
range of $p\left(x\right)$ over $I_{n}(\underline{d},\overline{d})$ using the
range enclosure property of Bernstein polynomials as follows:

###### Theorem 1 (Theorem 2 (Garloff, 1985a)).

Let $p$ be a multivariate polynomial of degree $L$ over the $n$-dimensional
box $I_{n}(\underline{d},\overline{d})$ with Bernstein coefficients $b_{K,L}$,
$0\leq K\leq L$. Then, for all $x\in I_{n}$, the following inequality holds:

(16) $\displaystyle\min_{K\leq L}\leavevmode\nobreak\ b_{K,L}\leq
p\left(x\right)\leq\max_{K\leq L}\leavevmode\nobreak\ b_{K,L}.$

The traditional approach to computing the range enclosure of $p$ is to compute
all the Bernstein coefficients of $p$ to determine their minimum and maximum
(Garloff, 1985b; Garloff and Smith, 2004; Zettler and Garloff, 1998). However,
computing all the coefficients has a complexity of
$\mathcal{O}\left(\left(l_{max}+1\right)^{n}\right)$, where
$l_{max}=\max\limits_{1\leq i\leq n}l_{i}$, which increases exponentially with
the dimension $n$. Luckily, the Bernstein coefficients enjoy monotonicity
properties, whenever the region $I_{n}(\underline{d},\overline{d})$ is
restricted to be an orthant (i.e., the sign of $x_{i}$ does not change within
$I_{n}(\underline{d},\overline{d})$, for each $i\in\\{1,\ldots,n\\}$) (Smith,
2009). Using such monotonicity properties, one can compute the minimum and
maximum Bernstein coefficients (denoted by $\underline{B}_{p,L}$,
$\overline{B}_{p,L}$) with a time complexity of
$\mathcal{O}\left(2\left(l_{max}+1\right)^{2}\right)$ which does not depend on
the dimension $n$.

###  5.2. Zero Crossing Estimation using only a few Bernstein Coefficients

Now we discuss how to use the range enclosure property above to reduce the
number of computed Bernstein coefficients. First, we note that the zero
crossing of a polynomial $p$ in a given input region $I_{n}$ depends on its
estimate range given by $\underline{B}_{p,L}$ and $\overline{B}_{p,L}$. More
specifically, if $\underline{B}_{p,L}>0$ ($\overline{B}_{p,L}<0$), then the
entire polynomial is positive (negative), which means that there is no zero-
crossing. If $\underline{B}_{p,L}$ and $\overline{B}_{p,L}$ have different
signs, and because of the estimation error of these bounds, the polynomial $p$
may still be positive, negative, or have a zero crossing in the region. In
this case, we need additional information such as the bounds of the gradient
of the polynomial $p$ within the input region, that are given by
$\underline{B}_{\nabla p,L}$ and $\overline{B}_{\nabla p,L}$ (which can be
computed efficiently thanks to the fact that gradients of polynomials are
polynomials themselves). Such additional information about the worst-case
gradient of the polynomial leads to a natural estimate of whether a zero
crossing occurs in a region.

Due to space constraints, we omit the analysis of bounding the estimation
error introduced by relying only on the maximum and minimum of the polynomial
$\underline{B}_{p,L}$ and $\overline{B}_{p,L}$ along with the maximum and
minimum of the gradient $\underline{B}_{\nabla p,L}$ and $\overline{B}_{\nabla
p,L}$. Instead, we support our claim using the empirical evidence shown in
Table 1. Using the same benchmarks used in Section 4.2, we train a third
neural network that takes as input only the four inputs
$\underline{B}_{p,L},\overline{B}_{p,L},\underline{B}_{\nabla
p,L},\overline{B}_{\nabla p,L}$ and compare its generalization performance
(column 7 of Table 1). As shown in the table, the third neural network
sacrifices some accuracy compared to the ones that use all Bernstein
coefficients. But on the other side, it reduces the overhead to compute the
Bernstein coefficients by order of magnitude as can be seen by comparing the
“execution overhead” reported in columns 4, 6, and 8 for the power basis, the
Bernstein basis, and the reduced Bernstein basis, respectively.

###  5.3. Search Space pruning using Bernstein Coefficients

The range enclosure property and the discussion above open the door for a
natural solution of the “alternative abstraction” problem mentioned in Section
3.2. The maximum and minimum Bernstein coefficients can be used as an
abstraction (in addition to convex upper and lower bounds) of high-order
polynomials. Such abstractions can be refined with every iteration of the
solver. They can be used to identify portions of the search space for which
one of the polynomials is guaranteed to be positive (and hence a solution does
not exist). More details about integrating this abstraction and the convex
abstraction are given in the implementation section below.

##  6\. Algorithm Architecture and Implementation Details

Algorithm 1 PolyARBerNN

Input: $I_{n}(\underline{d},\overline{d}),p_{1},p_{2},\ldots,p_{m},\epsilon$,
Output: $x_{\text{Sol}}$

1: $orthants:=\texttt{Partition\\_Region}(I_{n})$

2: $Neg:=\\{\\}$

3: $Ambig:=\\{orthants\\}$

4: $\text{List\\_pols}:=\\{p_{1},\ldots,p_{m}\\}$

5: while $\texttt{Compute\\_Maximum\\_Volume}(Ambig)\geq\epsilon$ do

6: $p:=\texttt{Select\\_Poly}\left(\text{List\\_pols},Neg\right)$

7:
$region:=\texttt{Remove\\_Ambiguous\\_Region\\_From\\_List}\left(Ambig\right)$

8: $\big{(}\underline{B}_{p,L},\overline{B}_{p,L},\underline{B}_{\nabla
p,L},\overline{B}_{\nabla p,L}\big{)}$ :=
$\texttt{Compute\\_Bern\\_Coeff}(p,region)$

9: if $\underline{B}_{p,L}>0$ then

10: break

11: else if $\underline{B}_{p,L}<0$ then

12: $Neg:=Neg\cup(p,L^{-}_{0}\left(p\right))$

13: break

14: end if

15: $\text{(under\\_approx, over\\_approx,
split)}:=NN(\underline{B}_{p,L},\overline{B}_{p,L},\underline{B}_{\nabla
p,L},\overline{B}_{\nabla p,L},region\big{)}$

16: $action:=\texttt{Select\\_best\\_action}(\text{under\\_approx,
over\\_approx, split})$

17:
$L^{-}_{0}\left(p\right),L^{+}_{0}\left(p\right),L^{+/-}_{0}\left(p\right):=\texttt{Convex\\_Abst\\_Refin\\_PolyAR}\left(p,action,region\right)$

18: $Ambig:=Amibg\cup L^{+/-}_{0}\left(p\right)$

19: $Neg:=Neg\cup(p,L^{-}_{0}\left(p\right))$

20: end while

21: if $\texttt{is\\_List\\_Empty}(Ambig)$ then

22: if A negative region in $Neg$ has all the polynomials then

23: $x_{\text{Sol}}:=$ any point in the negative region

24: else

25: return the problem is UNSAT

26: end if

27: else

28:
$x_{\text{Sol}}:=\texttt{CAD\\_Solver\\_Parallel}\left(Ambig,p_{1},\ldots,p_{m}\right)$

29: end if

In this section, we describe the implementation details of our solver
PolyARBerNN. As a pre-processing step, the tool divides the input region
$I_{n}$ into several regions such that each one is an orthant. This allows the
tool to process each orthant in parallel or sequentially. The tool keeps track
of all regions for which the sign of a polynomial is not fixed. These regions
are called ambiguous regions, and they are stored in a list called $Ambig$. As
long as the volume of the regions in this list is larger than a user-defined
threshold $\epsilon$, then our tool will continuously use abstractions to
identify portions in which one of the polynomials is always positive (and
hence removed from the search space) or negative (and hence the tool will give
higher priority for this region). The abstraction refinement is iteratively
applied in Lines 5-17 of Algorithm 1. In each abstraction refinement step, the
tool picks a polynomial $p$ and a region $region$ based on several heuristics
(Lines 6-7). In lines 8-14, we compute the maximum/minimum Bernstein
coefficients followed by checking the sign of the polynomial within this
region. Suppose the Bernstein coefficients indicate that the polynomial is
always positive in this region. In that case, this provides a guarantee that a
solution does not exist in this region (recall that Problem 1 searches for a
point where _all_ polynomials are negative). Similarly, if the polynomial is
always negative, then it will be added to the list of negative regions. For
those polynomials for which the Bernstein abstraction failed to identify their
signs, we query the trained neural network to estimate the best convex
abstraction possible (Lines 15-16). Based on the neural network suggestion, we
use the PolyAR tool (Yasser Shoukry, 2021) to compute the convex abstraction
(Line 17), which returns portions of this region that are guaranteed to belong
to the zero sublevel set $L^{-}_{0}(p)$, those who belong to the zero
superlevel set $L^{+}_{0}(p)$, and those remain ambiguous $L^{+/-}_{0}(p)$.
The process of using Bernstein abstraction and the convex abstraction (which
is guided by the trained neural network) continues until all remaining
ambiguous regions are smaller than a user-defined threshold $\epsilon$ in
which case it will be processed in parallel using a sound and complete tool
that implements Cylindrical Algebraic Decomposition (CAD) such as Z3 and Yices
(line 28 in Algorithm 1).

Figure 3. The architecture of the trained NN that is used to guide the
abstraction refinement process within PolyARBerNN. We used a fully connected
NN that contains an input layer with 4 neurons, three hidden layers with 40
neurons each, and one output layer with three neurons. All neurons are ReLU-
based except for the output neurons which uses SoftMax non-linearity.

The neural network itself is trained using randomly generated, quadratic, two-
dimensional polynomials where the coefficients follow a uniform distribution
between $-1$ and $1$. For each randomly generated polynomial, we used PolyAR
to compute the volumes of the $L^{+}_{0}(p),L^{-}_{0}(p),L^{+/-}_{0}(p)$
regions. We use a fully connected NN that contains an input layer, three
hidden layers, and one output layer (shown in Figure 3). The input layer has
four neurons, the hidden layers have 40 neurons each, and the output layer has
three neurons. We use a dropout of probability $0.5$ in the first and second
hidden layers to avoid overfitting. We use the ReLU activation function for
all the hidden layers and the Softmax activation function for the output
layer. We use Adam as an optimizer and cross-entropy as a loss function.
Although the neural network is trained on simple quadratic two-dimensional
polynomials, we observed it generalizes well to higher-order polynomials with
several variables. This will become apparent during the numerical evaluation
in which polynomials of different orders and several variables will be used to
evaluate the tool.

Correctness Guarantees: We conclude our discussion with the following result
which captures the correctness guarantees of the proposed tool:

###### Theorem 1.

The PolyARBerNN solver is sound and complete.

###### Proof.

This result follows from the fact that search space is pruned using sound
abstractions (convex upper bounds or Bernstein-based). The neural network and
the convex lower bound polynomials are just used as heuristics to guide the
refinement process. Finally, CAD-based algorithms (which are sound and
complete) are used to process the portions of the search space which are not
pruned by the abstraction refinement. ∎

##  7\. Generalization to polynomial optimization problems:

In this section, we focus on providing a solution to Problem 2. Our approach
is to turn the optimization problem (Problem 2) into a feasibility problem
(Problem 1). First, we recall that the gradient of $p$, $\nabla
p=[\frac{\partial p}{\partial x_{1}},\cdots{},\frac{\partial p}{\partial
x_{n}}]$, where $\frac{\partial p}{\partial x_{i}}$ is the partial derivative
of $p$ with respect to $x_{i}$, is a vector of $n$ polynomials. The optimal
value of $p$ occurs either (i) when the vector of partial derivatives are all
equal to zero or (ii) at the boundaries of the input region.

To find the critical points $x^{*}$ of $p$ where $\nabla
p\left(x^{*}\right)=0$, we add the $n$ polynomial constraints $\frac{\partial
p}{\partial x_{i}}\leq 0,1\leq i\leq n$ and $-\frac{\partial p}{\partial
x_{i}}\leq 0,1\leq i\leq n$ to the constraints of the optimization problem.
Now, we modify the PolyARBernNN solver to output _all_ possible regions in
which all the constraints are satisfied. This can be easily computed by taking
the intersections within the regions stored in the data structure $Neg$ in
Algorithm 1. These regions enjoy the property that _all_ points in these
regions are critical points of $p$. In addition, we modify PolyARBernNN to
output _all_ the remaining ambiguous regions whose volumes are smaller than
the user-specified threshold $\epsilon$ and for which the CAD-based solvers
returned a solution. These regions enjoy the property that _there exists_ a
point inside these regions which is a critical point. These modifications are
captured in (Line 2 of Algorithm 2).

Since the minimum/maximum of $p$ may occur at the boundaries of the region
$I_{n}\left(\underline{d},\overline{d}\right)$, our solver samples from the
boundaries of the region $I_{n}\left(\underline{d},\overline{d}\right)$ (Line
4 in Algorithm 2). The solver uses $2\sqrt{n}(\epsilon)^{1/n}$ as sampling
distance between two successive boundary samples—recall $\epsilon$ is a user-
defined parameter and was used in Algorithm 1 as a threshold on the refinement
process. Finally, we evaluate the polynomial $p$ in the obtained samples and
we take the minimum and maximum over the obtained values (Line 7 in Algorithm
2). All the details can be found in Algorithm 2.

Algorithm 2 PolyAROpt

Input: $I_{n}(\underline{d},\overline{d}),p,p_{1},p_{2},\ldots,p_{m},\epsilon$

Output: $\hat{p}_{min},\hat{p}_{max}$

1: $\nabla p=\texttt{Grad\\_Poly}(p)$

2:
$\hat{reg}_{\text{list}}=\texttt{PolyARBerNN}(I_{n}(\underline{d},\overline{d}),\nabla
p,p_{1},\ldots,p_{m},\epsilon)$

3: $\hat{x}_{\text{list}}=\texttt{center}(\hat{reg}_{\text{list}})$

4:
$x^{end}_{list}=\texttt{Sample\\_boundaries}(I_{n}(\underline{d},\overline{d}),\epsilon)$

5: $\hat{p}_{list}=p(\hat{x}_{\text{list}})$; $\quad
p^{end}_{list}=p(x^{end}_{\text{list}})$

6: $p_{list}=\hat{p}_{list}\cup p^{end}_{list}$

7: $\hat{x}_{min}=\arg\min(p_{list}$); $\quad\hat{x}_{max}=\arg\max(p_{list})$

8: $\hat{p}_{min}=p(\hat{x}_{min})$; $\quad\hat{p}_{max}=p(\hat{x}_{max})$

We conclude our discussion with the following result which captures the error
between the solutions provided by PolyAROpt and the global optima.

###### Theorem 1.

Let $p^{*}_{min}$ and $p^{*}_{max}$ be the global optimal points for the
solution of Problem 2. The solution obtained by Algorithm 2, denoted by
$\hat{p}_{min}$ and $\hat{p}_{max}$ satisfies the following:

(17) $\displaystyle\left\lVert\hat{p}_{min}-p^{*}_{min}\right\rVert\leq
2\omega_{p}\sqrt{n}(\epsilon)^{1/n},$
$\displaystyle\left\lVert\hat{p}_{max}-p^{*}_{max}\right\rVert\leq
2\omega_{p}\sqrt{n}(\epsilon)^{1/n}$

where $\omega_{p}$ is the Lipschitz constant of the polynomial $p$ and
$\epsilon>0$ is a user-defined error.

###### Proof.

We note that there are three cases that Algorithm 2 uses to compute the
critical points, $\hat{x}_{min}$ and $\hat{x}_{max}$, which corresponds to the
approximation of the minimum $\hat{p}_{min}$ and maximum $\hat{p}_{max}$ of
the polynomial:

  1. (1)

Using the center of the regions in the $Neg$ list

  2. (2)

Using the center of the regions in the $Ambig$ list

  3. (3)

Using samples from the boundaries

We proceed by case analysis. Case 1: First, we note that _all_ the points
within the $Neg$ regions satisfy that $\nabla p=0$ and hence the value of the
polynomial takes the same exact value overall the region, hence the value of
$p$ at the center $\hat{x}$ of the region is the same at the global optima
$x^{*}$. Case 2 and Case 3: If we can show that $\hat{x}_{min}$ and
$\hat{x}_{max}$ are bounded from the actual optimal points $x^{*}_{min}$ and
$x^{*}_{max}$ by:

(18) $\displaystyle\left\lVert\hat{x}_{min}-x^{*}_{min}\right\rVert\leq
2\sqrt{n}(\epsilon)^{1/n}$
$\displaystyle\left\lVert\hat{x}_{max}-x^{*}_{max}\right\rVert\leq
2\sqrt{n}(\epsilon)^{1/n}$

then (17) will follow directly from the definition of the Lipschitz continuity
of polynomials. To show that inequalities (18) hold we proceed by case
analysis. However (18) follow directly in Case 2 from the fact that $Ambig$
regions have a volume that is smaller than $\epsilon$ (Line 7 in Algorithm 1)
and hence the distance between any two points within the regions is bounded by
$2\sqrt{n}(\epsilon)^{1/n}$. Similarly, Case 3 follows from the fact that
Algorithm 2 samples from the boundaries with a maximum distance between the
samples that is equal to $2\sqrt{n}(\epsilon)^{1/n}$.

∎

##  8\. Numerical Results - NN Training

In this section, we show the details of training and evaluating the NN used to
help PolyARBerNN selecting the best convex abstraction. We evaluate the
trained NN on six different benchmarks. The benchmarks are different than the
training benchmarks with respect to the input region, the degree of the
polynomial, and the number of variables of the polynomial. All the experiments
were executed on an Intel Core i7 2.6-GHz processor with 16 GB of memory.

###  8.1. Training data collection and pre-processing

####  8.1.1. Data collection

To collect the data, we generated random quadratic two-dimensional
polynomials:

$q\left(x_{1},\leavevmode\nobreak\ x_{2}\right)\leavevmode\nobreak\
=\leavevmode\nobreak\
c_{1}x_{1}^{2}+c_{2}x_{2}+c_{3}x_{1}x_{2}+c_{4}x_{2}+c_{5}y+c_{6},$

where the coefficients $c_{1},\ldots,c_{6}$ follow a uniform distribution
between $-1$ and $1$.The random generated polynomials are defined over the
domain $I_{2}=\big{[}-2,2\big{]}^{2}$. For each randomly generated polynomial,
we perform the abstraction refinement on the domain $I_{2}$, iteratively. In
every iteration, we perform under-approximation, over-approximation of the
original polynomial over a selected ambiguous region, and a split of the
ambiguous region. Next, we compute the volume of the remaining ambiguous
region after each action was implemented. The labels are a one-hot vector of
dimension three where each component represents the action that leads to the
maximum reduction in the volume of the ambiguous region, either under-
approximation, over-approximation or divide the region into two regions. We
ran the abstraction refinement process on all the generated polynomials to
collect the data
$\big{(}\underline{B^{i}}_{p_{j},L},\overline{B^{i}}_{p_{j},L},\underline{B^{i}}_{\nabla
p_{j},L},\overline{B^{i}}_{\nabla p_{j},L}\big{)}$, where $i$ denote the index
of the sample. We generate $50000$ samples for training, $10000$ samples for
validation, and $10000$ for testing.

####  8.1.2. Data Normalization

In the literature of NN (LeCun et al., 2012), it is important to normalize the
data when the data vary across a wide range of values. This normalization
leads to faster training and improves the generalization performance of the NN
(LeCun et al., 2012). Therefore, we normalize all the input data to a zero
mean and unit variance by adopting a simple affine transformation
$\text{data\\_sample}\leftarrow\frac{\text{data\\_sample}-\mu}{\sigma}$, where
$\mu$ and $\sigma$ are the mean of the data and its standard deviation. The
$\mu$ and $\sigma$ parameters are initialized with respectively the empirical
mean and standard deviation of the dataset and they are computed offline
before the training.

Figure 4. Percentage in reduction of the volume of ambiguous regions along
with the NN output number (the number is at the top of histograms) for 20
samples for 6 evaluation benchmarks described in Table II.

###  8.2. NN’s evaluation

We evaluated the NN on 6 different benchmarks as follows:

  * •

In the first and second benchmarks, we generate the same random quadratic
polynomials but in a different domain: $\big{[}-4,4\big{]}^{2}$ and
$\big{[}-10,10\big{]}^{2}$. This choice is made to test the generalization of
the NN outside the data domains that were used in its training. This is
important since the Bernstein coefficients of a polynomial (the input to the
NN) depend on the input region $I_{n}$.

  * •

In the third and fourth benchmarks, we generated random polynomials with
degrees $4$ and $10$ over the domain $\big{[}-2,2\big{]}^{2}$. These
benchmarks are used to validate the generalization of the NN to polynomials of
orders higher than the ones used in its training.

  * •

Finally, in the fifth and sixth benchmarks, we generated random polynomials
with higher dimensions, i.e., with dimension $n=4$ and $n=7$.

In summary, these benchmarks will help us to answer the following question:
can the trained NN generalize to new data with different domains (benchmarks 1
and 2), higher orders (benchmarks 3 and 4), and higher dimensions (benchmarks
5 and 6)? More detail about the different benchmarks is shown in Table 2.

Figure 4 shows the performance of the trained NN over 20 random samples of
each of the six benchmarks. For each sample, we used the framework in (Yasser
Shoukry, 2021) to compute the ground-truth percentage in the reduction of the
volume of ambiguous regions after applying every action under-approximation,
over-approximation, or split. We then evaluated the NN on each sample and
reported in Figure 4 both the ground-truth reduction of the ambiguous regions
(as bars) against the index of the action suggested by the NN (as the text
above the bars). As it can be seen from Figure 4, except for the second sample
of the first benchmark, the NN outputs represent the actions that lead to the
maximum reduction of the ambiguous region’s volume.

Table 2. Evaluation of the trained NN on the six different benchmarks.

Benchmark $p\left(x\right)$ $n$ order region Accuracy 1
$c_{1}x_{1}^{2}+c_{2}x_{2}+c_{3}x_{1}x_{2}+c_{4}x_{2}+c_{5}y+c_{6}$ $2$ $2$
$[-4,4]^{2}$ $95\%$ 2
$c_{1}x_{1}^{2}+c_{2}x_{2}+c_{3}x_{1}x_{2}+c_{4}x_{2}+c_{5}y+c_{6}$ $2$ $2$
$[-10,10]^{2}$ $93\%$ 3
$c_{1}x_{1}^{4}+c_{2}x_{2}^{3}+c_{3}x_{1}^{4}x_{2}^{3}+c_{4}x_{1}^{3}+c_{5}$
$2$ $4$ $[-2,2]^{2}$ $88\%$ 4
$c_{1}x_{1}^{10}+c_{2}x_{2}^{5}+c_{3}x_{1}^{5}x_{2}^{3}+c_{4}x_{1}^{5}+c_{5}$
$2$ $10$ $[-2,2]^{2}$ $87\%$ 5
$c_{1}x_{1}^{3}+c_{2}x_{2}^{3}+c_{3}x_{3}^{3}+c_{4}x_{4}^{3}$ $4$ $3$
$[-2,2]^{4}$ $81\%$ 6
$c_{1}x_{1}^{3}+c_{2}x_{2}^{3}+c_{3}x_{3}^{3}+c_{4}x_{4}^{3}+c_{5}x_{5}^{3}+c_{6}x_{6}^{3}+c_{7}x_{7}^{3}$
$7$ $3$ $[-2,2]^{7}$ $80\%$

Finally, we ran the same experiment for 1000 samples and report the percentage
of samples for which the NN was able to predict the action that leads to the
maximum reduction in the ambiguous region’s volume. As it can be seen from
Table 2, the trained NN is able to generalize on the different benchmarks. For
instance, evaluating the NN on different domains results in the lowest
accuracy of $93\%$. Furthermore, evaluating the NN on polynomials with higher-
order results in an accuracy of $87\%$. Finally, the NN achieves $80\%$ on
higher dimension benchmarks.

$m$ | SAT/UNSAT | Execution times vs Polynomial Order | Execution times vs Number of Variables
---|---|---|---
|
1 | UNSAT |  |
SAT |  |
5 | UNSAT |  |
SAT |  |
10 | UNSAT |  |
SAT |  |

Figure 5. Scalability results of PolyARBerNN in the UNSAT case for $1$, $5$,
and $10$ constraints. (left) evolution of the execution time in seconds as a
function of the order of the polynomials, (right) evolution of the execution
time in seconds as a function of the number of variables. The timeout is equal
to 1 hour.

##  9\. Numerical Results - Scalability Results

In this section, we study the scalability of PolyARBerNN in terms of execution
times by varying the order, the number of variables, and the number of the
polynomial constraints for instances when a solution exists (the problem is
Satisfiabile or SAT for short) and when a solution does not exist (or UNSAT
for short). We will perform this study in comparison with state-of-the-art
solvers including our previous solver PolyAR, Z3 8.9, and Yices 2.6. Next, we
compare the performance of PolyARBerNN against a theorem prover named PVS
which implements a Bernstein library to solve multivariate polynomial
constraints (Munoz and Narkawicz, 2013). Finally, we compare the scalability
of the PolyAROpt optimizer against the built-in optimization library in Z3 8.9
to solve an unconstrained multivariate polynomial optimization problem with
varying order and number of variables.

###  9.1. Scalability of PolyARBerNN against other solvers

In this experiment, we compare the execution times of PolyARBerNN against the
PolyAR tool (Yasser Shoukry, 2021), Z3 8.9, and Yices 2.6. We consider two
instances of Problem 1: an UNSAT and SAT problems. For each instance, we
consider three scenarios, $m=1$, $m=5$, and $m=10$ where $m$ is the number of
polynomial constraints. First, we vary the order of the polynomials from 0 to
1000 while fixing the number of variables (and hence the dimension of the
search space) to two. Alternatively, we also fix the order of the polynomials
to $30$ while varying the number of variables from 1 to 200. We set the
timeout of the simulations to be $1$ hour. Figure 5 reports the execution
times for all the experiments whenever the problem is UNSAT and SAT. As
evidenced by the figures, PolyARBerNN succeeded to solve the instances of
Problem 1 for all orders and numbers of variables in a few seconds. For
instance, solving $10$ polynomial constraints with $200$ variables and a
maximum order of $30$ took around $20\leavevmode\nobreak\ s$ leading to a
speed-up of $200\times$ compared to Z3 and Yices. On the other hand, other
solvers are incapable of solving the polynomial constraints for all orders or
number of variables and they time out after one hour. These results show the
scalability of the proposed approach by including Bernstein coefficients to
prune the search space and a NN to guide the abstraction refinement.

###  9.2. Scalability of PolyARBerNN versus PVS

In this experiment, we want to investigate the effect of the NN on the overall
performance of the system. We compare PolyARBerNN against the PVS theorem
prover which also uses Bernstein coefficients to reason about polynomial
constraints. The polynomials and the domains of the associated variables are
given in (Munoz and Narkawicz, 2013) which were originally used to assess the
scalability of the PVS theorem prover. As evinced by the results in Table 3,
using a combination of Bernstein and adding the NN to guide the convex
abstractions (in PolyARBerNN) leads to additional savings in the execution
time, leading to a speedup of $483\times$ in the Heart dipole example.

Table 3. Scalability results for PolyARBerNN against PVS.

Problem PVS PolyARBerNN Times (sec) Speedup Schwefel $1.23$ $\mathbf{0.022}$
$54.9\times$ Reaction diffusion $0.25$ $\mathbf{0.018}$ $12.88\times$ Caprasse
$0.03$ $\mathbf{0.025}$ $0.19\times$ Lotka-Volterra $0.22$ $\mathbf{0.026}$
$7.46\times$ Butcher $0.46$ $\mathbf{0.029}$ $14.86\times$ Magnetism $1.82$
$\mathbf{0.028}$ $64.0\times$ Heart dipole $15.01$ $\mathbf{0.031}$
$483.19\times$

Figure 6. Scalability results of PolyAROpt for unconstrained optimization.
(left) evolution of the execution time in seconds as a function of the order
of the polynomial. (right) evolution of the execution time in seconds as a
function of the number of variables. The timeout is equal to 1 hour.

###  9.3. Scalability of PolyAROpt against other solvers

We compare the scalability results of PolyAROpt with the Z3 solver due to the
fact that Z3 has a built-in optimization library. Unfortunately, Yices does
not have such an optimizer. We set the timeout of the simulation to be 1 hour.
Figure 6 reports the execution times of two experiments of computing the
minimum and maximum of unconstrained optimization. As evidenced by the two
figures, PolyAROpt succeeded to solve the unconstrained optimization problem
for all orders and number of variables. For instance, solving the
unconstrained optimization problem with $70$ variables and a maximum order of
$3$ took around $50\leavevmode\nobreak\ seconds$. On the other hand, the Z3
solver is incapable of solving the unconstrained optimization problem for all
orders or number of variables and times out.

##  10\. Numerical Results - Use Cases

In this section, we provide two engineering use cases. The first one focuses
on the use of PolyARBerNN to synthesize stabilizing non-parametric controllers
for nonlinear dynamical systems. The second use case focuses on the use of
PolyAROPT to perform reachability analysis of polynomial dynamical systems.

###  10.1. Use Case 1: Nonlinear Controller Design for a Duffing Oscillator

In this subsection, we assess the scalability of the PolyARBerNN solver
compared to state-of-the-art solvers for synthesizing a non-parametric
controller for a Duffing oscillator reported by (Fotiou et al., 2006). All the
details of the dynamics of the oscillator and how we generated the polynomial
constraints can be found in (Yasser Shoukry, 2021). We denote by $n$ the
dimension of the Duffing oscillator. We consider two instances of the
controller synthesis problem for the Duffing oscillator with the following
parameters:

  * •

$n=3$, $\zeta=1.0$, $x\left(0\right)=[0.15,0.15,0.15]$,
$L_{1}\left(x\left(k\right),u\left(k\right)\right)=\left(-x_{1}^{3}\left(k\right)+x_{3}^{3}\left(k\right)+u\left(k\right)-2\right)^{51}$,
$L_{2}\left(x\left(k\right),u\left(k\right)\right)=x_{1}^{51}\left(k\right)x_{3}^{7}\left(k\right)+x_{1}^{9}\left(k\right)x_{3}^{5}\left(k\right)-5x_{2}^{4}\left(k\right)-x_{2}^{2}\left(k\right)u^{2}\left(k\right)$,
which results in $9$ polynomial constraints with $4$ variables and max
polynomial order of $153$.

  * •

$n=4$, $\zeta=1.75$, $x\left(0\right)=[0.1,0.1,0.01,0.1]$,
$L_{1}\left(x\left(k\right),u\left(k\right)\right)=x_{1}^{4}\left(k\right)+x_{2}^{4}\left(k\right)+x_{3}^{4}\left(k\right)+x_{4}^{4}\left(k\right)-u^{4}\left(k\right)$,
$L_{2}\left(x\left(k\right),u\left(k\right)\right)=-x_{1}^{51}\left(k\right)x_{3}^{20}\left(k\right)-5x_{2}^{4}\left(k\right)-x_{2}^{2}\left(k\right)u^{2}\left(k\right)$,
$L_{3}\left(x\left(k\right),u\left(k\right)\right)=\left(x_{1}x_{2}^{2}-u\left(k\right)-100\right)^{41}$,
which results in $12$ polynomial constraints with $5$ variables and max
polynomial order of $82$.

$n$ | State Space | Execution Time Evolution over time
---|---|---
3 |  |
4 |  |

Figure 7. Results of controlling the Duffing oscillator with different $n$
(left) evolution of the states $x_{1}(k)$ and $x_{2}(k)$ for the solvers in
the state-space, (right) evolution of the execution time of solvers during the
$12$ seconds. The timeout is equal to $60s$. Trajectories are truncated once
the solver exceeds the timeout limit.

We feed the resultant polynomial inequality constraint to PolyARBerNN, Yices,
and Z3. We solve the feasibility problem for $n=3$ and $n=4$. We set the
timeout to be $60s$. Figure 7 (left) shows the state-space evolution of the
controlled Duffing oscillator for different solvers for number of variables
$n$ of $3$ and $4$. Figure 7 (right) shows the evolution of the execution time
of the solvers during the $12$ seconds. As it can be seen from Figure 7, our
solver PolyARBerNN succeeded to find a control input $u$ that regulates the
state to the origin for all $n$. However, off-the-shelf solvers are incapable
of solving all two instances, and they early time out after
$60\leavevmode\nobreak\ seconds$ out of the simulated $12\leavevmode\nobreak\
seconds$. This shows the scalability of the proposed approach.

Figure 8. Reachability computation for the FitzHugh-Nagumo neuron model. Left:
using PolyAROpt for number of steps $K=50$. Right: using Sapo for number of
steps $K=50$.

Figure 9. Volume of the reachable set of states that are obtained using
PolyAROpt and Sapo (left) FitzHugh-Nagumo, (center) Duffing oscillator, and
(right) Jet flight.

###  10.2. Use Case 2: Reachability analysis of a discrete polynomial
dynamical systems

In this section, we show how to use PolyAROpt to compute the reachable set of
states for discrete-time polynomial systems. We consider a discrete polynomial
dynamical system of the following form:

(19) $\displaystyle x_{k+1}=f\left(x_{k}\right),\leavevmode\nobreak\
k\in\mathbb{N},x_{0}\in Q_{0},$

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a multivariate polynomial
map of a maximum degree $L=\left(l_{1},\cdots,l_{n}\right)$, and $Q_{0}$ is a
bounded polyhedron in $\mathbb{R}^{n}$. In this subsection, we consider a
bounded-time reachability analysis of the system in (19). Computing the exact
reachability sets of this type of dynamic system is hard. Therefore, we
overapproximate the exact set with a simplified set such as bounded polyhedra.
Bounded polyhedra are easy to handle and analyze. Computing the reachability
set, after a finite time $K$, involves computing sequentially the reachability
set at every time step $k$ using the following relation
$Q_{k+1}=f\left(Q_{k}\right),k=0,\cdots,K-1$, where $Q_{k}$ is a bounded
polyhedra. A bounded polyhedra is represented with an H-representation
$Q_{k}=\left(A_{k},b_{k}\right)=\\{x\in\mathbb{R}^{n}|A_{k}x\leq b_{k}\\}$,
where the inequality is a point-wise inequality. The template $A_{k}$
represents the directions of $Q_{k}$’s faces and $b_{k}$ represents their
positions. Given a polyhedra $Q_{k}=\left(A_{k},b_{k}\right)$, we need to
compute $Q_{k+1}=f\left(Q_{k}\right)$. We assume that
$A_{k}\in\mathbb{R}^{m\times n}$ is given. Now, we need to compute
$b_{k+1}\in\mathbb{R}^{m}$ which can be obtained through the following
optimization problem (Ben Sassi et al., 2012):

(20) $\displaystyle-b_{k+1,i}\leq u_{k+1,i}=\min\limits_{x\in
Q_{k}}-A_{k+1,i}f\left(x\right),\forall i=1,\cdots,m,$

where $A_{k+1,i}$ and $b_{k+1,i}$ are the $i^{th}$ row and component of the
templates $A_{k+1}$ and $b_{k+1}$, respectively. In every step
$k\in\mathbb{N}$, we compute an upper bound for $-b_{k+1,i}$, by solving the
optimization problem (20) using PolyAROpt and then an overapproximation of the
reachability set $Q_{k+1}$ is computed. In order to use PolyAROpt to solve
(20), we need to overapproximate the polyhedra $Q_{k}$ with a hyperreactangle
$R_{k}$. Therefore, the optimization problem (20) is modified as follows:

$\displaystyle-b_{k+1,i}\leq u_{k+1,i}=$ $\displaystyle\min\limits_{x\in
R_{k}}-A_{k+1,i}f\left(x\right),\forall i=1,\cdots,m,$ (21)
$\displaystyle\text{subject to}\leavevmode\nobreak\ x\leavevmode\nobreak\
\in\leavevmode\nobreak\ Q_{k}.$

We implemented the reachability computation method and tested it on three
dynamical systems:

  * •

FitzHugh-Nagumo Neuron: is a polynomial dynamic systems that models the
electrical activity of a neuron (FitzHugh, 1961). We performed reachability
analysis for $K=50$ time steps with the initial set of states
$Q_{0}=[0.9,1.1]\times[2.4,2.6]$,

  * •

Duffing Oscillator: is a discrete-time version of a nonlinear oscillator model
(Yasser Shoukry, 2021). We performed reachability analysis for $K=50$ time
steps with the initial set of states $Q_{0}=[2.49,2.51]\times[1.49,1.51]$,

  * •

Jet Flight: is a discrete-time version of a jet flight model (Chen, 2015). We
performed reachability analysis for $K=50$ time steps with the initial set of
states $Q_{0}=[0.9,1.2]\times[0.9,1.2]$.

All the details of the dynamics of the three dynamical systems and how we
generated the polynomial constraints can be found in (Ben Sassi et al., 2012;
FitzHugh, 1961; Yasser Shoukry, 2021; Chen, 2015). We computed the reachable
sets for these dynamical systems and compared our results with Sapo (Dreossi,
2017). Sapo is a tool proposed for the reachability analysis of the class of
discrete-time polynomial dynamical systems. Sapo linearizes the optimization
problem (20) using the Bernstein form of the polynomial and was shown to
outperform state-of-the-art reachability analysis tools like Flow∗ (Chen et
al., 2013). Figure 8 shows the reachable sets computed by PolyAROpt compared
to Sapo for the FitzHugh-Nagumo Neuron model. Inspecting the results in Figure
8 qualitatively shows that PolyAROpt is capable of computing tighter sets. To
quantitatively compare the results of PolyAROpt and Sapo, we compute the
volume of each reachable set (for different time steps). Figure 9 shows these
volumes for all the three dynamical systems mentioned above. As evident by the
results in Figure 9, PolyAROpt results in reachable sets that are tighter than
the one obtained from the Sapo thanks to PolyAROpt’s ability of solving the
polynomial optimization problem without any relaxations. Such ability to avoid
relaxation results in several orders of magnitude reduction in the volume of
the reachable sets compared to Sapo; a significant improvement in the analysis
of such dynamical systems.

Conclusions. In this paper, we proposed PolyARBerNN, a solver for polynomial
inequality constraints. We proposed a systematic methodology to design neural
networks that can be used to guide the abstraction refinement process by
bridging the gap between neural network properties and the properties of
polynomial representations. We showed that the use of Bernstein coefficients
leads the way to designing better neural network guides and provides an
additional abstraction that can be used to accelerate the solver. We
generalized the solver to reason about optimization problems. We demonstrated
that the proposed solver outperforms state-of-the-art tools by several orders
of magnitude.

###### Acknowledgements.

This work was supported by the National Science Foundation under grant numbers
#2002405 and #2139781.

## References

  * (1)
  * mor (2018) 2018\.  Under-approximating reach sets for polynomial continuous systems, author=Xue, Bai and Fränzle, Martin and Zhan, Naijun. In _Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control (part of CPS Week)_. 51–60.
  * An et al. (2018) Jie An, Naijun Zhan, Xiaoshan Li, Miaomiao Zhang, and Wang Yi. 2018.  Model checking bounded continuous-time extended linear duration invariants. In _Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control (part of CPS Week)_. 81–90.
  * Bak et al. (2019) Stanley Bak, Hoang-Dung Tran, and Taylor T Johnson. 2019\.  Numerical verification of affine systems with up to a billion dimensions. In _Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control_. 23–32.
  * Bello et al. (2016) Irwan Bello, Hieu Pham, Quoc V Le, Mohammad Norouzi, and Samy Bengio. 2016.  Neural combinatorial optimization with reinforcement learning.  _arXiv preprint arXiv:1611.09940_ (2016).
  * Ben Sassi et al. (2012) Mohamed Amin Ben Sassi, Romain Testylier, Thao Dang, and Antoine Girard. 2012.  Reachability analysis of polynomial systems using linear programming relaxations. In _International Symposium on Automated Technology for Verification and Analysis_. Springer, 137–151.
  * Boreale (2018) Michele Boreale. 2018\.  Algorithms for exact and approximate linear abstractions of polynomial continuous systems. In _Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control (part of CPS Week)_. 207–216.
  * Brisebarre et al. (2012) Nicolas Brisebarre, Mioara Joldeş, Érik Martin-Dorel, Micaela Mayero, Jean-Michel Muller, Ioana Paşca, Laurence Rideau, and Laurent Théry. 2012.  Rigorous polynomial approximation using Taylor models in Coq. In _NASA Formal Methods Symposium_. Springer, 85–99.
  * Brown (2001) Christopher W Brown. 2001\.  Improved projection for cylindrical algebraic decomposition.  _Journal of Symbolic Computation_ 32, 5 (2001), 447–465.
  * Chen (2015) Xin Chen. 2015.  _Reachability analysis of non-linear hybrid systems using taylor models_.  Ph. D. Dissertation. Fachgruppe Informatik, RWTH Aachen University.
  * Chen et al. (2013) Xin Chen, Erika Ábrahám, and Sriram Sankaranarayanan. 2013\.  Flow*: An analyzer for non-linear hybrid systems. In _International Conference on Computer Aided Verification_. Springer, 258–263.
  * Collins (1975) George E Collins. 1975\.  Quantifier elimination for real closed fields by cylindrical algebraic decomposition. In _Automata Theory and Formal Languages 2nd GI Conference Kaisers lautern, May 20–23, 1975_. Springer, 134–183.
  * De Moura and Bjørner (2008) Leonardo De Moura and Nikolaj Bjørner. 2008.  Z3: An efficient SMT solver. In _International conference on Tools and Algorithms for the Construction and Analysis of Systems_. 337–340.
  * Devlin et al. (2017) Jacob Devlin, Jonathan Uesato, Surya Bhupatiraju, Rishabh Singh, Abdel-rahman Mohamed, and Pushmeet Kohli. 2017.  Robustfill: Neural program learning under noisy i/o. In _International conference on machine learning_. PMLR, 990–998.
  * Dreossi (2017) Tommaso Dreossi. 2017\.  Sapo: Reachability computation and parameter synthesis of polynomial dynamical systems. In _Proceedings of the 20th International Conference on Hybrid Systems: Computation and Control_. 29–34.
  * Dutertre (2014) Bruno Dutertre. 2014\.  Yices 2.2. In _International Conference on Computer Aided Verification_. Springer, 737–744.
  * Dutta et al. (2019) Souradeep Dutta, Xin Chen, and Sriram Sankaranarayanan. 2019\.  Reachability analysis for neural feedback systems using regressive polynomial rule inference. In _Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control_. 157–168.
  * England and Davenport (2016) Matthew England and James H Davenport. 2016.  The complexity of cylindrical algebraic decomposition with respect to polynomial degree. In _International Workshop on Computer Algebra in Scientific Computing_. Springer, 172–192.
  * Farouki and Rajan (1987) Rida T Farouki and VT Rajan. 1987.  On the numerical condition of polynomials in Bernstein form.  _Computer Aided Geometric Design_ 4, 3 (1987), 191–216.
  * Ferlez and Shoukry (2020) James Ferlez and Yasser Shoukry. 2020.  AReN: assured ReLU NN architecture for model predictive control of LTI systems. In _Proceedings of the 23rd International Conference on Hybrid Systems: Computation and Control_. 1–11.
  * FitzHugh (1961) Richard FitzHugh. 1961\.  Impulses and physiological states in theoretical models of nerve membrane.  _Biophysical journal_ 1, 6 (1961), 445–466.
  * Fotiou et al. (2006) Ioannis A Fotiou, Philipp Rostalski, Pablo A Parrilo, and Manfred Morari. 2006.  Parametric optimization and optimal control using algebraic geometry methods.  _Internat. J. Control_ 79, 11 (2006), 1340–1358.
  * Garloff (1985a) Jürgen Garloff. 1985a.  Convergent bounds for the range of multivariate polynomials. In _International Symposium on Interval Mathematics_. Springer, 37–56.
  * Garloff (1985b) Jürgen Garloff. 1985b.  Convergent bounds for the range of multivariate polynomials. In _International Symposium on Interval Mathematics_. Springer, 37–56.
  * Garloff and Smith (2004) Jürgen Garloff and Andrew P Smith. 2004.  An improved method for the computation of affine lower bound functions for polynomials.  In _Frontiers in Global Optimization_. Springer, 135–144.
  * Gasse et al. (2019) Maxime Gasse, Didier Chételat, Nicola Ferroni, Laurent Charlin, and Andrea Lodi. 2019\.  Exact Combinatorial Optimization with Graph Convolutional Neural Networks. In _Advances in Neural Information Processing Systems 32_.
  * Goldfarb and Idnani (1983) Donald Goldfarb and Ashok Idnani. 1983.  A numerically stable dual method for solving strictly convex quadratic programs.  _Mathematical programming_ 27, 1 (1983), 1–33.
  * Hong (1990) H. Hong. 1990.  An Improvement of the Projection Operator in Cylindrical Algebraic Decomposition. In _Proceedings of the International Symposium on Symbolic and Algebraic Computation_ (Tokyo, Japan) _(ISSAC ’90)_. Association for Computing Machinery, New York, NY, USA, 261–264.  https://doi.org/10.1145/96877.96943
  * Irving et al. (2016) Geoffrey Irving, Christian Szegedy, Alexander A Alemi, Niklas Eén, François Chollet, and Josef Urban. 2016.  Deepmath-deep sequence models for premise selection.  _Advances in Neural Information Processing Systems_ 29 (2016), 2235–2243.
  * Kaliszyk et al. (2017) Cezary Kaliszyk, François Chollet, and Christian Szegedy. 2017.  Holstep: A machine learning dataset for higher-order logic theorem proving.  _arXiv preprint arXiv:1703.00426_ (2017).
  * Kochdumper and Althoff (2020) Niklas Kochdumper and Matthias Althoff. 2020.  Reachability analysis for hybrid systems with nonlinear guard sets. In _Proceedings of the 23rd International Conference on Hybrid Systems: Computation and Control_. 1–10.
  * LeCun et al. (2012) Yann A LeCun, Léon Bottou, Genevieve B Orr, and Klaus-Robert Müller. 2012.  Efficient backprop.  In _Neural networks: Tricks of the trade_. Springer, 9–48.
  * Lestingi et al. (2020) Livia Lestingi, Mehrnoosh Askarpour, Marcello M Bersani, and Matteo Rossi. 2020.  Formal verification of human-robot interaction in healthcare scenarios. In _International Conference on Software Engineering and Formal Methods_. Springer, 303–324.
  * Mahboubi (2006) Assia Mahboubi. 2006\.  Programming and certifying a CAD algorithm in the Coq system. In _Dagstuhl Seminar Proceedings_. Schloss Dagstuhl-Leibniz-Zentrum für Informatik.
  * Manual (1987) CPLEX User’s Manual. 1987\.  Ibm ilog cplex optimization studio.  _Version_ 12 (1987), 1987–2018.
  * Marcucci and Tedrake (2019) Tobia Marcucci and Russ Tedrake. 2019.  Mixed-integer formulations for optimal control of piecewise-affine systems. In _Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control_. 230–239.
  * McCallum (1998) Scott McCallum. 1998\.  An improved projection operation for cylindrical algebraic decomposition.  In _Quantifier Elimination and Cylindrical Algebraic Decomposition_. Springer, 242–268.
  * Munoz and Narkawicz (2013) César Munoz and Anthony Narkawicz. 2013.  Formalization of Bernstein polynomials and applications to global optimization.  _Journal of Automated Reasoning_ 51, 2 (2013), 151–196.
  * Munoz (2015) Cesar A Munoz. 2015\.  Formal Methods in Air Traffic Management: The Case of Unmanned Aircraft Systems. In _International Colloquium on Theoretical Aspects of Computing (ICTAC 2015)_.
  * Narkawicz and Munoz (2012) Anthony Narkawicz and César A Munoz. 2012.  Formal Verification of Conflict Detection Algorithms for Arbitrary Trajectories.  _Reliab. Comput._ 17 (2012), 209–237.
  * Nelson and Pecheur (2002) Stacy D Nelson and Charles Pecheur. 2002.  Formal verification for a next-generation space shuttle. In _International Workshop on Formal Approaches to Agent-Based Systems_. Springer, 53–67.
  * Rabi (2020) Maben Rabi. 2020\.  Piece-wise analytic trajectory computation for polytopic switching between stable affine systems. In _Proceedings of the 23rd International Conference on Hybrid Systems: Computation and Control_. 1–11.
  * Ray and Nataraj (2012) Shashwati Ray and PSV Nataraj. 2012.  A Matrix Method for Efficient Computation of Bernstein Coefficients.  _Reliab. Comput._ 17, 1 (2012), 40–71.
  * Roy et al. (2018) Debayan Roy, Michael Balszun, Thomas Heurung, and Samarjit Chakraborty. 2018.  Multi-Domain Coupling for Automated Synthesis of Distributed Cyber-Physical Systems. In _2018 IEEE International Symposium on Circuits and Systems (ISCAS)_. IEEE, 1–5.
  * Sadraddini and Tedrake (2020) Sadra Sadraddini and Russ Tedrake. 2020.  Robust output feedback control with guaranteed constraint satisfaction. In _Proceedings of the 23rd International Conference on Hybrid Systems: Computation and Control_. 1–10.
  * Selsam et al. (2019) Daniel Selsam, Matthew Lamm, Benedikt Bünz, Percy Liang, Leonardo de Moura, and David L. Dill. 2019\.  Learning a SAT Solver from Single-Bit Supervision. In _International Conference on Learning Representations_.  https://openreview.net/forum?id=HJMC_iA5tm
  * Shen et al. (2019) Zuowei Shen, Haizhao Yang, and Shijun Zhang. 2019\.  Deep network approximation characterized by number of neurons.  _arXiv preprint arXiv:1906.05497_ (2019).
  * Shoukry et al. (2018) Yasser Shoukry, Michelle Chong, Masashi Wakaiki, Pierluigi Nuzzo, Alberto Sangiovanni-Vincentelli, Sanjit A Seshia, Joao P Hespanha, and Paulo Tabuada. 2018.  SMT-based observer design for cyber-physical systems under sensor attacks.  _ACM Transactions on Cyber-Physical Systems_ 2, 1 (2018), 1–27.
  * Smith (2009) Andrew Paul Smith. 2009\.  Fast construction of constant bound functions for sparse polynomials.  _Journal of Global Optimization_ 43, 2 (2009), 445–458.
  * Sun et al. (2019) Xiaowu Sun, Haitham Khedr, and Yasser Shoukry. 2019\.  Formal verification of neural network controlled autonomous systems. In _Proceedings of the 22nd ACM International Conference on Hybrid Systems: Computation and Control_. 147–156.
  * Vandenberghe (2010) Lieven Vandenberghe. 2010\.  The CVXOPT linear and quadratic cone program solvers.  _Online: http://cvxopt. org/documentation/coneprog. pdf_ (2010).
  * Vinod and Oishi (2018) Abraham P Vinod and Meeko MK Oishi. 2018.  Scalable underapproximative verification of stochastic LTI systems using convexity and compactness. In _Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control (Part of CPS Week)_. 1–10.
  * Xue et al. (2019) Bai Xue, Qiuye Wang, Naijun Zhan, and Martin Fränzle. 2019\.  Robust invariant sets generation for state-constrained perturbed polynomial systems. In _Proceedings of the 22nd ACM international conference on hybrid systems: Computation and control_. 128–137.
  * Yasser Shoukry (2021) Wael Fatnassi Yasser Shoukry. 2021\.  PolyAR: A Highly Parallelizable Solver For Polynomial Inequality Constraints Using Convex Abstraction Refinement. In _IFAC-PapersOnLine_ , Vol. 54. 43–48.  Issue 5.
  * Zettler and Garloff (1998) Michael Zettler and Jürgen Garloff. 1998.  Robustness analysis of polynomials with polynomial parameter dependency using Bernstein expansion.  _IEEE Trans. Automat. Control_ 43, 3 (1998), 425–431.

##  Appendix A Proof of Proposition 4

###### Proof.

Note that $NN_{a_{p}\rightarrow X_{0}}(a_{p})$ is a vector-valued function
which returns the roots $X_{0}=(x^{1}_{0},x^{2}_{0},\cdots,x^{n_{r}}_{0})$ of
the polynomial $p$. Therefore, to upper bound the Lipschitz constant of
$NN_{a_{p}\rightarrow X_{0}}(a_{p})$, we will start by upper bounding the
Lipschitz constant of its components functions $NN^{j}_{a_{p}\rightarrow
x^{j}_{0}}(a_{p})$, $1\leq j\leq n_{r}$. To that end, consider two polynomials
with coefficients $a_{p}$ and $a^{\prime}_{p}$ such that $\left\lVert
a_{p}-a^{\prime}_{p}\right\rVert\leq\epsilon$. Therefore:

(22) $\displaystyle\left\lVert NN^{j}_{a_{p}\rightarrow
x^{j}_{0}}(a_{p})-NN^{j}_{a_{p}\rightarrow
x^{j}_{0}}(a^{\prime}_{p})\right\rVert_{2}=\left\lVert
x^{j}_{0}(a_{p})-x^{j}_{0}(a^{\prime}_{p})\right\rVert_{2}\leq
C_{a_{p}}(x^{j}_{0})\left\lVert
a_{p}-a^{\prime}_{p}\right\rVert\leq\overline{C}_{a_{p}}\epsilon$

where $x^{j}_{0}(a_{p})$ and $x^{j}_{0}(a^{\prime}_{p})$ are the location of
the $j$th root for the polynomials with coefficients $a_{p}$ and
$a^{\prime}_{p}$, respectively. The last two inequalities follow from the
definition of the condition number (Definition 2). Now,

(23) $\displaystyle\left\lVert NN_{a_{p}\rightarrow
X_{0}}(a_{p})-NN_{a_{p}\rightarrow X_{0}}(a^{\prime}_{p})\right\rVert_{2}$
$\displaystyle=\sqrt{\sum\limits_{j=1}^{n_{r}}\left\lVert
NN^{j}_{a_{p}\rightarrow x^{j}_{0}}(a_{p})-NN^{j}_{a_{p}\rightarrow
x^{j}_{0}}(a^{\prime}_{p})\right\rVert_{2}^{2}}$ (24)
$\displaystyle\leq\sqrt{n_{r}\overline{C}_{a_{p}}^{2}\epsilon^{2}}=\sqrt{n_{r}}\overline{C}_{a_{p}}\epsilon.$

From which we conclude that the Lipschitz constant of $NN_{a_{p}\rightarrow
X_{0}}(a_{p})$ is bounded by
$\sqrt{n_{r}}\overline{C}_{a_{p}}=\mathcal{O}(n_{r}\overline{C}_{a_{p}})$. ∎

##  Appendix B Proof of Proposition 5

###### Proof.

We assume that the number of sub-regions $l$ is fixed for each dimension $n$,
and $l^{\frac{1}{n}}=k\leavevmode\nobreak\ \in\mathbb{N}$. Partitioning the
input space into $l$ sub-regions occurs by dividing the interval for each
dimension into $k$ sub-intervals. Without loss of generality, the sub-neural
network $NN_{I_{n}\rightarrow I^{i}_{n}}(I_{n})$ for $n=1$ can be defined as:

(25) $(\underline{d}^{i},\overline{d}^{i})=NN_{I_{n}\rightarrow
I^{i}_{n}}(I_{n})=NN_{I_{n}\rightarrow
I^{i}_{n}}(\underline{d},\overline{d})=\left(\underline{d}+\frac{i}{k}\left(\overline{d}-\underline{d}\right),\;\overline{d}+\frac{i+1}{k}\left(\overline{d}-\underline{d}\right)\right)=\begin{bmatrix}1+i/k&-i/k\\\
-(i+1)/k&1+(i+1)/k\end{bmatrix}\begin{bmatrix}\underline{d}\\\
\overline{d}\end{bmatrix}$

A generalization to a higher dimension is straightforward by replacing $i$
with a multi-index in each dimension. Note that $NN_{I_{n}\rightarrow
I^{i}_{n}}(I_{n})$ is a multivariate linear function in its inputs and hence
its Lipschitz constant can be computed as the largest singular value. Indeed,
the linear function depends only on the constant $k$ (which depends on the
constant $l$ and the dimension $n$) from which we conclude that the Lipschitz
constant of $NN_{I_{n}\rightarrow I^{i}_{n}}$ is some constant
$\mathcal{O}(n)$.

∎

##  Appendix C Proof of Proposition 6

###### Proof.

It follows from equations (2)-(4) that the sub-neural network
$NN_{X_{0}\rightarrow ZC_{i}}$ can be written as:

(26) $\displaystyle ZC_{i}(a_{p},I^{i}_{n})$
$\displaystyle=NN_{X_{0}\rightarrow ZC_{i}}\left(NN_{a_{p}\rightarrow
X_{0}}(a_{p}),NN_{I_{n}\rightarrow
I^{i}_{n}}(I_{n})\right)=NN_{X_{0}\rightarrow
ZC_{i}}\left((x_{0}^{1},\ldots,x_{0}^{n_{r}}),(\underline{d}^{i},\overline{d}^{i})\right).$

The indicator variable $ZC_{i}$ should be set to zero whenever all the roots
$x_{0}^{j}$ lies outside the hyperrectangle
$I_{n}^{i}(\underline{d}^{i},\overline{d}^{i})$. First note that a root
$x_{0}^{j}$ lies outside $I_{n}^{i}(\underline{d}^{i},\overline{d}^{i})$ if
and only if the following condition holds:

(27) $\displaystyle x^{j}_{0}\not\in
I_{n}^{i}(\underline{d}^{i},\overline{d}^{i})\quad\Longleftrightarrow\quad\sum\limits_{k=1}^{k=n}\bigg{|}x^{j}_{0,k}-\underline{d}^{i}_{k}\bigg{|}+\bigg{|}x^{j}_{0,k}-\overline{d}^{i}_{k}\bigg{|}-\sum\limits_{k=1}^{n}\left(\overline{d}^{i}_{k}-\underline{d}^{i}_{k}\right)>0$

where $x^{j}_{0,k},\underline{d}^{i}_{k},$ and $\overline{d}^{i}_{k}$ are the
$k$th element in the vectors $x^{j}_{0},\underline{d}^{i}$ and
$\overline{d}^{i}$, respectively. Hence, the indicator variable $ZC_{i}$
should be set to zero whenever the following conditions hold:

(28) $\displaystyle
ZC_{i}(a_{p},I^{i}_{n})=0\quad\Longleftrightarrow\quad\text{max}_{j\in\\{1,\ldots,n_{r}\\}}\left(\sum\limits_{k=1}^{k=n}\bigg{|}x^{j}_{0,k}-\underline{d}^{i}_{k}\bigg{|}+\bigg{|}x^{j}_{0,k}-\overline{d}^{i}_{k}\bigg{|}-\sum\limits_{k=1}^{n}\left(\overline{d}^{i}_{k}-\underline{d}^{i}_{k}\right)\right)=0$

Before we compute the Lipschitz constant of $NN_{X_{0}\rightarrow ZC_{i}}$ in
Equation (28), we recall the following identities. Consider two functions
$f(x)$ and $g(x)$ with Lipschitz constants $L_{f}$ and $L_{g}$, respectively.
Then:

  * •

The Lipschitz constant of $\max(f(x),g(x))$ is bounded by $L_{f}+L_{g}$.

  * •

The Lipschitz constant of $f(x)+g(x)$ is bounded by $\max(L_{f},L_{g})$.

Now notice that
$|x^{j}_{0,k}-\underline{d}^{i}_{k}|=max(x^{j}_{0,k}-\underline{d}^{i}_{k},0)+max(-x^{j}_{0,k}+\underline{d}^{i}_{k},0)$.
Applying the identities above along with the fact that the Lipschitz constant
of $x^{j}_{0,k}-\underline{d}^{i}_{k}$ is $\mathcal{O}(1)$, we conclude that
the Lipschitz constant of $|x^{j}_{0,k}-\underline{d}^{i}_{k}|$ is
$\mathcal{O}(1)$. Hence, the Lipschitz constant of
$\sum\limits_{k=1}^{k=n}\bigg{|}x^{j}_{0,k}-\underline{d}^{i}_{k}\bigg{|}+\bigg{|}x^{j}_{0,k}-\overline{d}^{i}_{k}\bigg{|}-\sum\limits_{k=1}^{n}\left(\overline{d}^{i}_{k}-\underline{d}^{i}_{k}\right)$
is also $\mathcal{O}(1)$. Finally, the Lipschitz constant of the right hand
side of Equation (28) is $\mathcal{O}(n_{r})$. We conclude our proof by
noticing that all the operators in Equation (28)—namely the absolute value,
the max operator, summation, and checking the final value against a
constant—can be implemented exactly using ReLU neural networks (Ferlez and
Shoukry, 2020) and hence the neural network $NN_{X_{0}\rightarrow ZC_{i}}$
will also have a Lipschitz constant equal to $\mathcal{O}(n_{r})$.

∎

##  Appendix D Proof of Proposition 7

###### Proof.

This result follows directly by noticing that $NN_{ZC\rightarrow L^{+}/L^{-}}$
can be computed as a linear function:

(29) $\displaystyle NN_{ZC\rightarrow
L^{+}/L^{-}}(ZC_{1},\ldots,ZC_{l})=v\sum_{i=1}^{l}(1-ZC_{i})$

where $v$ is a constant that depends on the volume of the hyperrecatngle
$I_{n}$. Since $l$ is a constant, we conclude the result.

∎