AMP short = AMP , long = asymmetric multicore processor BVP short = BVP , long
= boundary value problem CLEVER short = Clever , long = Cross Lipschitz
Extreme Value for nEtwork Robustness CPPO short = cppo , long = complete
pointed partial order CPS short = CPS , long = cyber-physical system DAG short
= DAG , long = directed acyclic graph DBF short = DBF , long = Digital Beam
Former DCPO short = dcpo , long = directed-complete partial order DSL short =
DSL , long = domain-specific language EDSL short = EDSL , long = embedded
domain-specific language FBP short = FBP , long = free boundary problem FivePs
short = 5Ps , long = five Ps HS short = HS , long = hybrid system ISA short =
ISA , long = instruction set architecture IVP short = IVP , long = initial
value problem ODE short = ODE , long = ordinary differential equation PC short
= PC , long = Pulse Compression PCROP short = PC-ROP , long = PDE-constrained
rearrangement optimization problem PDCPO short = pointed dcpo , long = pointed
directed-complete partial order PDE short = PDE , long = partial differential
equation PI short = PI , long = principal investigator POSET short = poset ,
long = partially ordered set RE short = r. e. , long = recursively enumerable
RNN short = RNN , long = recurrent neural network ROP short = ROP , long =
rearrangement optimization problem SAC short = Sac , long = Single Assignment
C SBSE short = SBSE , long = search-based software engineerign SOP short = SOP
, long = shape optimization problem TTE short = TTE , long = Type-II Theory of
Effectivity

# A Domain-Theoretic Framework for Robustness Analysis of Neural Networks

Can Zhou Department of Computer Science, University of Oxford, Oxford, United
Kingdom<EMAIL_ADDRESS>Razin A. Shaikh Department of Computer Science,
University of Oxford, Oxford, United Kingdom, and Quantinuum Ltd., Oxford,
United Kingdom<EMAIL_ADDRESS>Yiran Li School of Computer Science,
University of Nottingham Ningbo China<EMAIL_ADDRESS>Amin
Farjudian School of Computer Science, University of Nottingham Ningbo China,
<EMAIL_ADDRESS>(Corresponding author)

###### Abstract

A domain-theoretic framework is presented for validated robustness analysis of
neural networks. First, global robustness of a general class of networks is
analyzed. Then, using the fact that Edalat’s domain-theoretic $L$-derivative
coincides with Clarke’s generalized gradient, the framework is extended for
attack-agnostic local robustness analysis. The proposed framework is ideal for
designing algorithms which are correct by construction. This claim is
exemplified by developing a validated algorithm for estimation of Lipschitz
constant of feedforward regressors. The completeness of the algorithm is
proved over differentiable networks, and also over general position
$\operatorname{\mathrm{ReLU}}$ networks. Computability results are obtained
within the framework of effectively given domains. Using the proposed domain
model, differentiable and non-differentiable networks can be analyzed
uniformly. The validated algorithm is implemented using arbitrary-precision
interval arithmetic, and the results of some experiments are presented. The
software implementation is truly validated, as it handles floating-point
errors as well.

Keywords: domain theory, neural network, robustness, Lipschitz constant,
Clarke-gradient

MSC classes: 06B35, 68Q55, 49J52, 68T37

##  1 Introduction

A system is said to be robust with respect to perturbations of a set of
parameters, if small changes to those parameters do not result in significant
changes in the behavior of the system. Robustness is a core requirement of
safety-critical systems, including neural networks deployed in high-stakes
applications [18, 44]. As a result, several methods have been suggested for
_measuring_ robustness, e. g., through adversarial attacks [40, 9] or via
attack-agnostic methods [45, 85, 52, 49, 50]. A common and effective attack-
agnostic approach involves estimation of the local Lipschitz constant of the
network. This makes sense as larger Lipschitz constants signify more
sensitivity to input values. As such, a useful estimate of the Lipschitz
constant should be a _tight upper bound_. A lower bound of the Lipschitz
constant results in false negatives—which are more serious than false
positives—and a very loose upper bound generates too many false positives,
which renders the estimate ineffective. Sound and accurate computation of the
Lipschitz constant is also essential in Lipschitz regularization of neural
networks [4, 66].

Of particular interest is the estimation of Lipschitz constant of non-
differentiable networks, such as $\operatorname{\mathrm{ReLU}}$ networks.
Computing the Lipschitz constant of $\operatorname{\mathrm{ReLU}}$ networks is
not just an NP-hard problem [82], it is strongly inapproximable [49, Theorem
4]. As a result, methods of Lipschitz analysis fall into two _disjoint_
categories:

  1. (1)

Scalable, but not validated, e. g., [85, 12, 56];

  2. (2)

Validated, but not scalable, e. g., [45, 49, 50], and the method of the
current article.

The validated methods—though not scalable—are essential in providing the
ground truth for verification of the scalable methods. For instance, consider
the scalable CLEVER method of [85]. [49] demonstrate that, on certain parts of
the input domain, CLEVER outputs inaccurate results, even on small networks.
The method of [49, 50], in turn, although based on a validated theory, is
implemented using floating-point arithmetic, which is prone to inaccuracies,
as will be demonstrated in Section 5. It is well-known that floating-point
errors lead to vulnerability in safety-critical neural networks [48, 89].

Based on the theory of continuous domains, we develop a framework for
robustness analysis of neural networks, within which we design an algorithm
for estimation of local Lipschitz constant of feedforward regressors. We prove
the soundness of the algorithm. Moreover, we prove completeness of the
algorithm over differentiable networks, and also over general position
$\operatorname{\mathrm{ReLU}}$ networks.

To avoid the pitfalls of floating-point errors, we implement our algorithm
using arbitrary-precision interval arithmetic (with outward rounding) which
guarantees soundness of the implementation as well. As such, our method is
_fully validated_ , that is, not only does it follow a validated theory, but
also it has a validated implementation.

###  1.1 Related work

A common approach to robustness analysis of neural networks is based on
estimation of Lipschitz constants. A mathematical justification of this
reduction is provided in [45, Theorem 2.1] and [85, Theorem 3.2]. The global
Lipschitz constant of a network may be estimated by multiplying the Lipschitz
constants of individual layers. This approach, taken by [80], is scalable but
provides loose bounds. In CLEVER [85], a statistical approach is adopted using
Extreme Value Theory [42]. While CLEVER is scalable, it may underestimate the
Lipschitz constant significantly [49]. A sound approach to Lipschitz
estimation was introduced by [45], but the method applies only to
differentiable shallow networks. Semidefinite programming has been used by
[36] and [43] to obtain tight upper bounds on the Lipschitz constant of a
network.

Of particular relevance to our work are methods that are based on abstract
interpretation [15]. The study of abstract interpretation in a domain-
theoretic framework goes back (at least) to [8, 2]. In the current context,
the main idea behind abstract interpretation is to enclose values, and then
follow how these enclosures propagate through the layers of a given network.
Various types of enclosures may be used, e. g., hyperboxes, zonotopes, or
polytopes, in increasing order of accuracy, and decreasing order of efficiency
[37, 59, 74]. Theoretically, methods that are based on interval analysis [83,
67, 58] or convex outer polytopes [87, 76] are also subsumed by the framework
of abstract interpretation. These methods, however, are applied over the
network itself, rather than its gradient.

In contrast, in [11, 10], abstract interpretations were used for estimating
Lipschitz constants of programs with control flow, which include neural
networks as special cases. [49, 50] have used Clarke’s generalized gradient
[14] to analyze non-differentiable (e. g., $\operatorname{\mathrm{ReLU}}$)
networks. Their overapproximation of the Clarke-gradient also follows the
spirit of abstract interpretation. Based on dual numbers, sound automatic
differentiation, and Clarke Jacobian, [57] have improved upon the results of
[11] and [49, 50].

Domains and neural networks serve different purposes, and there has been
little interaction between the two. One notable exception is the use of
probabilistic power domains for analysis of forgetful Hopfield networks [26].
[55, 54] also present a domain-theoretic framework for analysis of a general
class of networks which communicate real values. Their focus is on certain
distributed systems which include neural networks as a special case, but they
do not discuss the issue of robustness. Domain theory, however, has been used
in robustness analysis, but not in the context of neural networks.
Specifically, [60] have developed a framework for robustness analysis, with an
emphasis on hybrid systems. Their work provides an important part of the
foundation of the current paper. The results of [60] were developed further in
[61, 62, 35], in the context of abstract robustness analysis.

As we will analyze maximization of interval functions in our domain-theoretic
framework, we also point out the functional algorithm developed by [73] for
maximization. The focus of [73] is the functional style of the algorithm,
although he uses a domain semantic model. The domain model used in [73] is,
however, algebraic, which is not suitable for our framework. Our models are
all non-algebraic continuous domains.

###  1.2 Contributions

At a theoretical level, the main contribution of the current work is a domain-
theoretic framework, which provides a unified foundation for validated methods
of robustness and security analysis of neural networks, including interval
method of [83], Clarke-gradient methods of [49, 50], [7], and [57], and those
based on abstract interpretation, e. g., [37, 74]. In terms of results, the
main contributions are as follows:

Global robustness analysis:

We prove that feedforward classifiers always have tightest robust
_approximations_ (Theorem 3.5), and regressors always have robust domain-
theoretic _extensions_ (Theorem 3.9).

Local robustness analysis:

Using Edalat’s domain-theoretic $L$-derivative—which coincides with Clarke-
gradient—we develop an algorithm for computing the Lipschitz constant of
feedforward networks. In particular, we prove _soundness and completeness_ of
the algorithm. We also prove _computability_ of our method within the
framework of effectively given domains. These results are presented in Lemma
4.12 and Theorems 4.13, 4.16, 4.19 and 4.21. To the best of our knowledge,
completeness and computability results are missing in the previous literature
on the subject.

Fully validated method:

As the domains involved are effectively given, our method has a simple and
direct implementation using arbitrary-precision interval arithmetic, which is
correct by construction. Furthermore, following the completeness theorems, we
can obtain the results to within any given degree of accuracy, for a broad
class of networks.

###  1.3 Structure of the paper

The rest of this article is structured as follows: The preliminaries and
technical background are presented in Section 2. Section 3 contains the
results on global robustness analysis. In Section 4, we present the main
results on validated local robustness analysis. Section 5 contains the account
of some relevant experiments. Section 6 contains the concluding remarks.

##  2 Preliminaries

In this section, we present the technical background which will be used later
in the paper.

###  2.1 Metric Structure

We must first specify the metric used in measuring perturbations of the input,
and the norm used on the relevant gradients. For every $n\in\mathbb{N}$ and
$p\in[1,\infty]$, let $\ell_{p}^{n}$ denote the normed space
$(\mathbb{R}^{n},{\left\lVert\thinspace{.}\thinspace\right\rVert}^{(n)}_{p})$,
in which the norm ${\left\lVert\thinspace{.}\thinspace\right\rVert}^{(n)}_{p}$
is defined for all $x=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$ by
${\left\lVert\thinspace{x}\thinspace\right\rVert}^{(n)}_{p}\coloneqq\left(\sum_{i=1}^{n}\lvert\thinspace{x_{i}}\thinspace\rvert^{p}\right)^{1/p}$,
if $p\in[1,\infty)$, and
${\left\lVert\thinspace{x}\thinspace\right\rVert}^{(n)}_{p}\coloneqq\max\\{{\lvert\thinspace{x_{i}}\thinspace\rvert}\mid{1\leq
i\leq n}\\}$, if $p=\infty$. We let
$d_{p}^{(n)}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to[0,\infty)$ denote the
metric induced by
${\left\lVert\thinspace{.}\thinspace\right\rVert}^{(n)}_{p}$, i. e.,
$d_{p}^{(n)}(x,y)={\left\lVert\thinspace{x-y}\thinspace\right\rVert}^{(n)}_{p}$.
When $n$ is clear from the context, we drop the superscript and write
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{p}$ and $d_{p}$. For every
$p\in[1,\infty]$, define the conjugate $p^{\prime}\in[1,\infty]$ as follows:

$p^{\prime}\coloneqq\left\\{\begin{array}[]{ll}\infty,&\text{if }p=1,\\\
p/(p-1),&\text{if }p\in(1,\infty),\\\ 1,&\text{if }p=\infty.\\\
\end{array}\right.$ (1)

In our discussion, whenever the perturbation of the input is measured using
${\left\lVert\thinspace{.}\thinspace\right\rVert}^{(n)}_{p}$-norm, then the
norm of the gradient must be taken to be
${\left\lVert\thinspace{.}\thinspace\right\rVert}^{(n)}_{p^{\prime}}$.
Furthermore, as the spaces $\ell^{n}_{p}$ are finite-dimensional, the topology
on these spaces does not depend on $p$, a fact which is not true in infinite-
dimensional $\ell_{p}$ spaces. Thus, in the current article, whenever we
mention a closed, open, or compact subset of $\mathbb{R}^{n}$, the topology is
assumed to be that of $\ell^{n}_{\infty}$.

######  Remark 2.1.

The readers familiar with functional analysis may notice that the Banach space
$\ell^{n}_{p^{\prime}}$ is the continuous dual of the Banach space
$\ell^{n}_{p}$. This explains the correspondence between the norm on the
gradient and the perturbation norm. To keep the discussion simple, we do not
present any technical background from functional analysis. The interested
reader may refer to classical texts on the subject, such as [70, 3]

######  Remark 2.2.

We will only consider perturbations with respect to the
${\left\lVert\thinspace{.}\thinspace\right\rVert}^{(n)}_{p}$-norm, which is
the default norm considered in most of the related literature. We point out,
however, that other metrics have also been suggested for quantifying
perturbations, most notably, the Wasserstein metric [88]. In fact, the choice
of an adequate metric remains an open problem [72].

###  2.2 Compact Subsets

Assume that $f:\mathbb{R}^{n}\to\mathbb{R}$ is (the semantic model of) a
regressor neural network with $n$ input neurons. In robustness analysis, we
must consider the action of $f$ over suitable _subsets_ of the input space.
Specifically, we work with closed subsets of both the input and output spaces.
Let us justify this choice formally.

Let $(\mathbb{S},d)$ be a metric space with distance
$d:\mathbb{S}\times\mathbb{S}\to\mathbb{R}$. For any $X\subseteq\mathbb{S}$
and $\delta>0$, we define the open neighborhood:

$B(X,\delta)\coloneqq\\{y\in\mathbb{S}\mid\exists x\in X:d(x,y)<\delta\\}.$
(2)

The set $B(X,\delta)$ is indeed open, because it is the union of open balls
$B(\left\\{{s}\right\\},\delta)$ with $s\in X$. Intuitively, $B(X,\delta)$ is
the set of points in $X$ with imprecision less than $\delta$. For all
$X\subseteq\mathbb{S}$ and $\delta>0$, we have
$X\subseteq\overline{X}\subseteq B(X,\delta)=B(\overline{X},\delta)$, in which
$\overline{X}$ denotes the closure of $X$. Thus, _in the presence of
imprecision, $X$ and $\overline{X}$ are indistinguishable_. Furthermore, it is
straightforward to verify that
$\overline{X}=\bigcap_{\delta>0}B(X,\delta)=\bigcap_{\delta>0}\overline{B(X,\delta)}$.
Intuitively, $\overline{B(X,\delta)}$ is the set of points in $X$ with
imprecision less than or equal to $\delta$. Thus, _the closure $\overline{X}$
is the set of points that are in $X$ with arbitrarily small imprecision_.

In practice, the input space of a network is a bounded region such as
$\mathbb{S}=[-M,M]^{n}$. As $[-M,M]^{n}$ is compact, and closed subsets of
compact spaces are compact, in this paper, we will be working with compact
subsets. A more detailed discussion of the choice of closed subsets, in the
broader context of robustness analysis, may be found in [62].

###  2.3 Domain Theory

Domain theory has its roots in topological algebra [39, 51], and it has
enriched computer science with powerful methods from order theory, algebra,
and topology. Domains gained prominence when they were introduced as a
mathematical model for lambda calculus by [71]. We have already mentioned that
our domain-theoretic framework provides a unified foundation for validated
methods of robustness and security analysis of neural networks. Domain theory,
however, provides more than just a unifying framework:

  1. (1)

Domain-theoretic models are ideal for designing algorithms which are _correct
by construction_ , as will be demonstrated in Section 4.

  2. (2)

Domains provide a denotational semantics for computable analysis according to
TTE [84, Theorem 9.5.2]. Our computability results (Theorems 4.16, 4.19, and
4.21) are obtained within the framework of effectively given domains [75]. As
a consequence, our constructs can be implemented directly on digital
computers.

We present a brief reminder of the concepts and notations that will be needed
later. The interested reader may refer to [39, 1, 38, 41] for more on domains
in general, and refer to [33, 27] for the interval domain, in particular. A
succinct and informative survey of the theory may be found in [51].

For any subset $X$ of a POSET $(D,\sqsubseteq)$, we write $\bigsqcup X$ and
$\bigsqcap X$ to denote the least upper bound, and the greatest lower bound,
of $X$, respectively, whenever they exist. In our discussion, $x\sqsubseteq y$
may be interpreted as ‘ _$y$ contains more information than $x$_’. A subset
$X\subseteq D$ is said to be _directed_ if $X\neq\emptyset$ and every two
elements of $X$ have an upper bound in $X$, i. e., $\forall x,y\in X:\exists
z\in X:(x\sqsubseteq z)\wedge(y\sqsubseteq z)$. An $\omega$-chain—i. e., a set
of the form $x_{0}\sqsubseteq x_{1}\sqsubseteq\cdots\sqsubseteq
x_{n}\sqsubseteq\cdots$—is an important special case of a directed set. A
POSET $(D,\sqsubseteq)$ is said to be a _DCPO_ if every directed subset
$X\subseteq D$ has a least upper bound in $D$. The POSET $(D,\sqsubseteq)$ is
said to be _pointed_ if it has a bottom element, which we usually denote by
$\bot_{D}$, or if $D$ is clear from the context, simply by $\bot$. Directed-
completeness is a basic requirement for almost all the POSET in our
discussion.

######  Example 2.3.

Assume that $\Sigma\coloneqq\left\\{{0,1}\right\\}$ and let $\Sigma^{*}$
denote the set of finite strings over $\Sigma$, ordered under the prefix
relation. Then, $\Sigma^{*}$ is a POSET, but it is not directed-complete. For
instance, the set $\\{{0^{n}}\mid{n\in\mathbb{N}}\\}$ is an $\omega$-chain,
but it does not have the least upper bound, or any upper bound. If we let
$\Sigma^{\omega}$ denote the set of countably infinite sequences over
$\Sigma$, then the POSET
$\Sigma^{\infty}\coloneqq\Sigma^{*}\cup\Sigma^{\omega}$, ordered under prefix
relation, is a DCPO, and may be regarded as a completion of $\Sigma^{*}$. The
POSET $\Sigma^{\infty}$ is indeed a PDCPO with the empty string as the bottom
element.

We let $\mathbb{K}\mathbb{R}^{n}_{\bot}$ denote the POSET with the carrier set
$\\{{C\subseteq\mathbb{R}^{n}}\mid{C\text{ is non-empty and
compact}}\\}\cup\left\\{{\mathbb{R}^{n}}\right\\}$, ordered by reverse
inclusion, i. e., $\forall X,Y\in\mathbb{K}\mathbb{R}^{n}_{\bot}:X\sqsubseteq
Y\iff Y\subseteq X$. By further requiring the subsets to be convex, we obtain
the sub-poset $\mathbb{C}\mathbb{R}^{n}_{\bot}$. Finally, we let
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$ denote the POSET of hyperboxes of
$\mathbb{R}^{n}$—i. e., subsets of the form
$\prod_{i=1}^{n}[a_{i},b_{i}]$—ordered under reverse inclusion, with
$\mathbb{R}^{n}$ added as the bottom element. The three POSET are PDCPO and
$\bigsqcup X=\bigcap X$, for any directed subset $X$.

######  Remark 2.4.

The reverse inclusion order $X\sqsubseteq Y\iff Y\subseteq X$ conforms to the
interpretation of $X\sqsubseteq Y$ as ‘$Y$ contains more information than
$X$’, in that, compared with the larger set $X$, the subset $Y$ is a more
accurate approximation of the result.

A central concept in domain theory is that of the _way-below_ relation. Assume
that $(D,\sqsubseteq)$ is a DCPO and let $x,y\in D$. The element $x$ is said
to be _way-below_ $y$—written as $x\ll y$—iff for every directed subset $X$ of
$D$, if $y\sqsubseteq\bigsqcup X$, then there exists an element $d\in X$ such
that $x\sqsubseteq d$. Informally, $x$ is way-below $y$ iff it is impossible
to get past $y$ (by suprema of directed sets) without first getting past $x$.
We may also interpret $x\ll y$ as ‘$x$ is a finitary approximation of $y$’, or
even as, ‘$x$ is a lot simpler than $y$’ [1, Section 2.2.1].

It is easy to prove that $\forall x,y\in D:x\ll y\implies x\sqsubseteq y$. The
converse, however, is not always true. For instance, over
$\mathbb{K}\mathbb{R}^{n}_{\bot}$ (hence, also over
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$ and $\mathbb{C}\mathbb{R}^{n}_{\bot}$) we
have the following characterization:

$\forall K_{1},K_{2}\in\mathbb{K}\mathbb{R}^{n}_{\bot}:K_{1}\ll K_{2}\iff
K_{2}\subseteq{K_{1}}^{\circ},$ (3)

in which ${K_{1}}^{\circ}$ denotes the interior of $K_{1}$.

######  Example 2.5.

Consider the interval domain ${\mathbb{I}}\mathbb{R}^{1}_{\bot}$ and let
$x\coloneqq[0,2]$ and $y\coloneqq[0,1]$. Although $x\sqsubseteq y$, we have
${x}^{\circ}=(0,2)$. Hence, by (3), we have $x\not\ll y$. We can also prove
that $x\not\ll y$ without referring to (3). Take the $\omega$-chain
$y_{n}\coloneqq[-2^{-n},1+2^{-n}]$ for $n\in\mathbb{N}$. We have
$y=\bigsqcup_{n\in\mathbb{N}}y_{n}$, but $\forall n\in\mathbb{N}:x\nsqsubseteq
y_{n}$.

######  Example 2.6.

In the poset $\Sigma^{\infty}$ of Example 2.3, we have: $x\ll
y\iff(x\sqsubseteq y)\wedge(x\in\Sigma^{*})$.

For every element $x$ of a DCPO $(D,\sqsubseteq)$, let
$\twoheaddownarrow{x}\coloneqq\\{{a\in D}\mid{a\ll x}\\}$. In domain-theoretic
terms, the elements of $\twoheaddownarrow{x}$ are the true approximants of
$x$. In fact, the way-below relation is also known as the order of
approximation [1, Section 2.2.1]. A subset $B$ of a DCPO $(D,\sqsubseteq)$ is
said to be a _basis_ for $D$, if for every element $x\in D$, the set
$B_{x}\coloneqq\twoheaddownarrow{x}\cap B$ is a directed subset with supremum
$x$, i. e., $x=\bigsqcup B_{x}$. A DCPO $(D,\sqsubseteq)$ is said to be
continuous if it has a basis, and $\omega$-continuous if it has a countable
basis. We call $(D,\sqsubseteq)$ a _domain_ if it is a continuous PDCPO. The
PDCPO ${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, $\mathbb{K}\mathbb{R}^{n}_{\bot}$,
and $\mathbb{C}\mathbb{R}^{n}_{\bot}$, are all $\omega$-continuous domains,
because:

  * •

$B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}\coloneqq\left\\{{\mathbb{R}^{n}}\right\\}\cup\\{{C\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}}\mid{C\text{
is a hyperbox with rational coordinates}}\\}$ is a basis for
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$;

  * •

$B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}\coloneqq\left\\{{\mathbb{R}^{n}}\right\\}\cup\\{{C\in\mathbb{K}\mathbb{R}^{n}_{\bot}}\mid{C\text{
is a finite union of hyperboxes with rational coordinates}}\\}$ is a basis for
$\mathbb{K}\mathbb{R}^{n}_{\bot}$;

  * •

$B_{\mathbb{C}\mathbb{R}^{n}_{\bot}}\coloneqq\left\\{{\mathbb{R}^{n}}\right\\}\cup\\{{C\in\mathbb{C}\mathbb{R}^{n}_{\bot}}\mid{C\text{
is a convex polytope with rational coordinates}}\\}$ is a basis for
$\mathbb{C}\mathbb{R}^{n}_{\bot}$.

The poset $\Sigma^{\infty}$ of Example 2.3 is also an $\omega$-continuous
domain, with $\Sigma^{*}$ as a countable basis. There is a significant
difference between $\Sigma^{\infty}$ and the domains
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, $\mathbb{K}\mathbb{R}^{n}_{\bot}$, and
$\mathbb{C}\mathbb{R}^{n}_{\bot}$, in that, for the elements of the basis
$\Sigma^{*}$ of $\Sigma^{\infty}$, we have: $\forall x\in\Sigma^{*}:x\ll x$.
By (3), this is clearly not true of the bases
$B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}$, $B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$
and $B_{\mathbb{C}\mathbb{R}^{n}_{\bot}}$, or indeed, any other bases for
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, $\mathbb{K}\mathbb{R}^{n}_{\bot}$, and
$\mathbb{C}\mathbb{R}^{n}_{\bot}$. As such, the domain $\Sigma^{\infty}$ is an
example of an algebraic domain, whereas the domains that we use in the current
paper are non-algebraic. Although algebraic domains have been used in real
number computation [16, 34] and, more broadly, representation of topological
spaces [79, 78], non-algebraic domains have proven more suitable for
computation over continuous spaces, e. g., in dynamical systems [20], exact
real number computation [33, 27], differential equation solving [22, 29], and
reachability analysis of hybrid systems [31, 60], to name a few. Hence, in
this article, we will be working in the framework of non-algebraic
$\omega$-continuous domains.

Apart from the order-theoretic structure, domains also have a topological
structure. Assume that $(D,\sqsubseteq)$ is a POSET. A subset $O\subseteq D$
is said to be _Scott open_ if it has the following properties:

  1. (1)

It is an upper set, that is, $\forall x\in O,\forall y\in D:x\sqsubseteq
y\implies y\in O$.

  2. (2)

It is inaccessible by suprema of directed sets, that is, for every directed
set $X\subseteq D$ for which $\bigsqcup X$ exists, if $\bigsqcup X\in O$ then
$X\cap O\neq\emptyset$.

######  Example 2.7.

Consider the poset $\mathbb{R}$ under the usual order on real numbers. Then,
for any $a\in\mathbb{R}$, the set $\\{{x\in\mathbb{R}}\mid{x>a}\\}$ is Scott
open. Note that $(\mathbb{R},\leq)$ is not a domain, because it is not even
complete. If we extend $\mathbb{R}$ to
$\mathbb{R}_{\pm\infty}\coloneqq[-\infty,+\infty]$ with $\forall
x\in\mathbb{R}:-\infty<x<+\infty$, then $(\mathbb{R}_{\pm\infty},\leq)$ is an
$\omega$-continuous domain (in fact, an $\omega$-continuous lattice) with
$\mathbb{Q}\cup\left\\{{-\infty}\right\\}$ as a basis, and each set of the
form $(a,+\infty]$ is Scott open.

######  Example 2.8.

Consider the domain ${\mathbb{I}}\mathbb{R}_{\bot}$. The set
$\\{{[a,b]}\mid{[a,b]\subseteq(0,1)}\\}$ is Scott open in
${\mathbb{I}}\mathbb{R}_{\bot}$. See also Proposition 2.10
LABEL:[below]prop:Scott_wayaboves.

The collection of all Scott open subsets of a POSET forms a $T_{0}$ topology,
referred to as the Scott topology. A function
$f:(D_{1},\sqsubseteq_{1})\to(D_{2},\sqsubseteq_{2})$ is said to be Scott-
continuous if it is continuous with respect to the Scott topologies on $D_{1}$
and $D_{2}$. Scott continuity can be stated purely in order-theoretic terms. A
map $f:(D_{1},\sqsubseteq_{1})\to(D_{2},\sqsubseteq_{2})$ between two posets
is Scott continuous if and only if it is monotonic and preserves the suprema
of directed sets, i. e., for every directed set $X\subseteq D_{1}$ for which
$\bigsqcup X$ exists, we have $f(\bigsqcup X)=\bigsqcup f(X)$ [41, Proposition
4.3.5]. When $D_{1}$ and $D_{2}$ are $\omega$-continuous domains, then we have
an even simpler formulation:

######  Proposition 2.9 ([1, Proposition 2.2.14]).

Assume that $f:D_{1}\to D_{2}$ is a map between two $\omega$-continuous
domains. Then, $f$ is Scott continuous iff it is monotonic and preserves the
suprema of $\omega$-chains.

For every element $x$ of a DCPO $(D,\sqsubseteq)$, let
$\twoheaduparrow{x}\coloneqq\\{{a\in D}\mid{x\ll a}\\}$.

######  Proposition 2.10 ([1, Proposition 2.3.6]).

Let $D$ be a domain with a basis $B$. Then, for each $x\in D$, the set
$\twoheaduparrow{x}$ is Scott open, and the collection
${\mathcal{O}}\coloneqq\left\\{{\twoheaduparrow{x}}\;\vrule\;{x\in
B}\right\\}$ forms a base for the Scott topology.

The maximal elements of ${\mathbb{I}}\mathbb{R}^{n}_{\bot}$,
$\mathbb{K}\mathbb{R}^{n}_{\bot}$, and $\mathbb{C}\mathbb{R}^{n}_{\bot}$ are
singletons, and the sets of maximal elements may be identified with
$\mathbb{R}^{n}$. For simplicity, we write $x$ to denote a maximal element
$\left\\{{x}\right\\}$. As a corollary of Proposition 2.10 and
characterization (3), we obtain:

######  Corollary 2.11.

Let ${\cal O}_{S}$ be the Scott topology on
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, $\mathbb{K}\mathbb{R}^{n}_{\bot}$, or
$\mathbb{C}\mathbb{R}^{n}_{\bot}$. Then, the restriction of
${\mathcal{O}}_{S}$ over $\mathbb{R}^{n}$ is the Euclidean topology.

Thus, the sets of maximal elements are homeomorphic to $\mathbb{R}^{n}$. For
any $K\geq 0$, by restricting to $[-K,K]^{n}$, we obtain the
$\omega$-continuous domains ${\mathbb{I}}[-K,K]^{n}$, $\mathbb{K}[-K,K]^{n}$,
and $\mathbb{C}[-K,K]^{n}$, respectively.

We say that a POSET $(D,\sqsubseteq)$ is _bounded-complete_ if each bounded
pair $x,y\in D$ has a supremum. Assume that $(X,\Omega(X))$ is a topological
space, and let $(D,\sqsubseteq_{D})$ be a bounded-complete domain. We let
$(D,\Sigma({D}))$ denote the topological space with carrier set $D$ under the
Scott topology $\Sigma({D})$. The space $[X\to D]$ of functions $f:X\to D$
which are $(\Omega(X),\Sigma({D}))$ continuous can be ordered pointwise by
defining:

$\forall f,g\in[X\to D]:\quad f\sqsubseteq g\iff\forall x\in
X:f(x)\sqsubseteq_{D}g(x).$

It is straightfoward to verify that the POSET $([X\to D],\sqsubseteq)$ is
directed-complete and $\forall x\in X:(\bigsqcup_{i\in
I}f_{i})(x)=\bigsqcup\\{{f_{i}(x)}\mid{i\in I}\\}$, for any directed subset
$\\{{f_{i}}\mid{i\in I}\\}$ of $[X\to D]$. For the POSET $([X\to
D],\sqsubseteq)$ to be continuous, however, the topological space
$(X,\Omega(X))$ must be _core-compact_ , as we explain briefly. Consider the
POSET $(\Omega(X),\subseteq)$ of open subsets of $X$ ordered under subset
relation. For any topological space $X$, this POSET is a complete lattice,
with $\emptyset$ as the bottom element, and $X$ as the top element.
Furthermore, we have:

$\forall A\subseteq\Omega(X):\quad\bigsqcup A=\bigcup A\text{ and }\bigsqcap
A={(\bigcap A)}^{\circ}.$

A topological space $(X,\Omega(X))$ is said to be core-compact if and only if
the lattice $(\Omega(X),\subseteq)$ is continuous.

######  Theorem 2.12.

For any topological space $(X,\Omega(X))$ and non-singleton bounded-complete
continuous domain $(D,\sqsubseteq_{D})$, the function space $([X\to
D],\sqsubseteq)$ is a bounded-complete continuous domain $\iff(X,\Omega(X))$
is core-compact.

###### Proof.

For the ($\Leftarrow$) direction, see [32, Proposition 2]. A proof of the
($\Rightarrow$) direction can also be found on [32, pages 62 and 63]. ∎

######  Notation 2.13 ($X\Rightarrow D$).

For any core-compact topological space $(X,\Omega(X))$ and bounded-complete
continuous domain $(D,\sqsubseteq_{D})$, we denote the continuous domain
$([X\to D],\sqsubseteq)$ by the notation $X\Rightarrow D$.

All the domains that will be used in the framework developed in this
article—including, ${\mathbb{I}}[-K,K]^{n}$, $\mathbb{K}[-K,K]^{n}$, and
$\mathbb{C}[-K,K]^{n}$—are bounded-complete, and they are core-compact under
their respective Scott topologies. Furthermore, the space $[-K,K]^{n}$ under
the Euclidean topology is core-compact. For more on core-compact spaces, the
interested reader may refer to, e. g., [39, Chapter II] and [41, Chapter 5].

Throughout this article, by an _interval function_ over a set $X$ we mean a
function $f:X\to{\mathbb{I}}\mathbb{R}^{n}_{\bot}$ for some $n\geq 1$. For any
set $X\subseteq\mathbb{R}^{n}$, we let $D^{(0)}(X)$ denote the continuous
domain of interval functions $X\Rightarrow{\mathbb{I}}\mathbb{R}_{\bot}$. Each
function in $D^{(0)}(X)$ is Euclidean-Scott-continuous, i. e., continuous with
respect to the Euclidean topology on $X$ and Scott topology on
${\mathbb{I}}\mathbb{R}_{\bot}$. The functions in $D^{(0)}(X)$ have a useful
characterization in terms of semicontinuity. Recall that, if $X$ is any
topological space, then:

  * •

$f:X\to\mathbb{R}$ is said to be upper semi-continuous at $x_{0}\in X$ iff,
for every $y>f(x_{0})$, there exists a neighborhood $U$ of $x_{0}$ (in the
topology of $X$) such that $\forall x\in U:f(x)<y$.

  * •

$f:X\to\mathbb{R}$ is said to be lower semi-continuous at $x_{0}\in X$ iff,
for every $y<f(x_{0})$, there exists a neighborhood $U$ of $x_{0}$ such that
$\forall x\in U:f(x)>y$.

  * •

$f:X\to\mathbb{R}$ is said to be upper (respectively, lower) semi-continuous
iff it is upper (respectively, lower) semi-continuous at every $x_{0}\in X$.

######  Proposition 2.14 ([30]).

Assume that $X\subseteq\mathbb{R}^{n}$ and
$f\equiv[{\underline{f}},{\overline{f}}]\in D^{(0)}(X)$, in which
${\underline{f}}$ and ${\overline{f}}$ are the lower and upper bounds of the
interval function $f$, respectively. Then, $f$ is Euclidean-Scott-continuous
$\iff$ ${\underline{f}}$ is lower semicontinuous, and ${\overline{f}}$ is
upper semicontinuous.

For any compact set $K\in\mathbb{K}\mathbb{R}^{n}_{\bot}$, we define the
hyperbox closure as:

$K^{\Box}\coloneqq\bigsqcup{\\{{R\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}}\mid{K\subseteq
R}\\}}.$ (4)

As $B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}$ is a basis for
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, it is straightforward to show that:

######  Proposition 2.15.

For any compact set $K\in\mathbb{K}\mathbb{R}^{n}_{\bot}$, we have
$K^{\Box}=\bigsqcup{\\{{b\in B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}}\mid{b\ll
K}\\}}$.

For each $i\in\left\\{{1,\ldots,n}\right\\}$, we let
$\pi_{i}:\mathbb{R}^{n}\to\mathbb{R}$ denote projection over the $i$-th
component, that is, $\pi_{i}(x_{1},\ldots,x_{n})\coloneqq x_{i}$. The
projections are all continuous functions, so they map compact sets to compact
sets. Hence, for any compact set $K\in\mathbb{K}\mathbb{R}^{n}_{\bot}$, the
set $\pi_{i}(K)$ is a compact subset of $\mathbb{R}$, for any given
$i\in\left\\{{1,\ldots,n}\right\\}$. For each
$i\in\left\\{{1,\ldots,n}\right\\}$, we define
${\underline{K_{i}}}\coloneqq\min\pi_{i}(K)$ and
${\overline{K_{i}}}\coloneqq\max\pi_{i}(K)$. It is straightforward to verify
that:

$\forall K\in\mathbb{K}\mathbb{R}^{n}_{\bot}:\quad
K^{\Box}=\prod_{i=1}^{n}[{\underline{K_{i}}},{\overline{K_{i}}}].$ (5)

######  Lemma 2.16.

For any hyperbox $b\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}$ and compact set
$K\in\mathbb{K}\mathbb{R}^{n}_{\bot}$, we have: $b\ll K\iff b\ll K^{\Box}$.

###### Proof.

Clearly $b\ll K^{\Box}$ implies $b\ll K$. For the opposite direction, assume
that $b=\prod_{i=1}^{n}b_{i}$. From (3), we deduce that $b\ll K$ iff $\forall
i\in\left\\{{1,\ldots,n}\right\\}:({\underline{b_{i}}}<{\underline{K_{i}}})\wedge({\overline{K_{i}}}<{\overline{b_{i}}})$.
This, combined with (5), entails $b\ll K^{\Box}$. ∎

######  Corollary 2.17.

The map
$(\cdot)^{\Box}:\mathbb{K}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{n}_{\bot}$
is Scott-continuous.

###### Proof.

As both $\mathbb{K}\mathbb{R}^{n}_{\bot}$ and
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$ are $\omega$-continuous, by Proposition
2.9, it suffices to show that the map $(\cdot)^{\Box}$ is monotonic and
preserves the suprema of $\omega$-chains. Monotonicity is straightforward.
Next, assume that $(K_{i})_{i\in\mathbb{N}}$ is an $\omega$-chain in
$\mathbb{K}\mathbb{R}^{n}_{\bot}$. We must show that
$\left(\bigsqcup_{i\in\mathbb{N}}K_{i}\right)^{\Box}=\bigsqcup_{i\in\mathbb{N}}K_{i}^{\Box}$.
The $\sqsupseteq$ direction follows from monotonicity. To prove the
$\sqsubseteq$ direction, by Proposition 2.15, we have:

$(\bigsqcup_{i\in\mathbb{N}}K_{i})^{\Box}=\bigsqcup{\left\\{{b\in
B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}}\;\vrule\;{b\ll\bigsqcup_{i\in\mathbb{N}}K_{i}}\right\\}}.$

Take any arbitrary $b\in B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}$ satisfying
$b\ll\bigsqcup_{i\in\mathbb{N}}K_{i}$. By Proposition 2.10, the set
$\twoheaduparrow{b}$ is Scott-open, which entails that $\exists
i_{0}\in\mathbb{N}:b\ll K_{i_{0}}$. By Lemma 2.16, we must have $b\ll
K_{i_{0}}^{\Box}$, from which we deduce that
$b\ll\bigsqcup_{i\in\mathbb{N}}K_{i}^{\Box}$. As $b$ was chosen arbitrarily,
then we must have
$\left(\bigsqcup_{i\in\mathbb{N}}K_{i}\right)^{\Box}\sqsubseteq\bigsqcup_{i\in\mathbb{N}}K_{i}^{\Box}$.
∎

######  Corollary 2.18.

For every $n\geq 1$, ${\mathbb{I}}\mathbb{R}^{n}_{\bot}$ is a sub-domain of
$\mathbb{K}\mathbb{R}^{n}_{\bot}$.

###### Proof.

Consider the map
$\iota:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to\mathbb{K}\mathbb{R}^{n}_{\bot}$
defined by $\forall x\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}:\iota(x)\coloneqq
x$. It is straightforward to prove that $\iota$ is Scott-continuous,
$\iota\circ(\cdot)^{\Box}\sqsubseteq\text{id}_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$,
and $(\cdot)^{\Box}\circ\iota=\text{id}_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}$.
Therefore, $(\iota,(\cdot)^{\Box})$ form a continuous embedding-projection
pair. ∎

######  Definition 2.19 (Extension, Canonical Interval Extension
${\boldsymbol{I}}f$, Approximation).

  1. (i)

A map $u:\mathbb{K}\mathbb{R}^{n}_{\bot}\to\mathbb{K}\mathbb{R}^{m}_{\bot}$ is
said to be an extension of $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ iff $\forall
x\in\mathbb{R}^{n}:u(\left\\{{x}\right\\})=\left\\{{f(x)}\right\\}$.

  2. (ii)

A map
$u:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{m}_{\bot}$ is
said to be an interval extension of
$f:\mathbb{R}^{n}\to\mathbb{K}\mathbb{R}^{m}_{\bot}$ iff $\forall
x\in\mathbb{R}^{n}:u(\left\\{{x}\right\\})=f(x)^{\Box}$.

  3. (iii)

For any $f:\mathbb{R}^{n}\to\mathbb{K}\mathbb{R}^{m}_{\bot}$, we define the
canonical interval extension
${{\boldsymbol{I}}f}:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{m}_{\bot}$
by:

$\forall\alpha\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}:\quad{{\boldsymbol{I}}f}(\alpha)\coloneqq\bigsqcap_{x\in\alpha}f(x)^{\Box}.$
(6)

  4. (iv)

A map
$u:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{m}_{\bot}$ is
said to be an interval approximation of
$f:\mathbb{R}^{n}\to\mathbb{K}\mathbb{R}^{m}_{\bot}$ if
$u\sqsubseteq{\boldsymbol{I}}f$.

######  Proposition 2.20.

For every Euclidean-Scott-continuous
$f:\mathbb{R}^{n}\to\mathbb{K}\mathbb{R}^{m}_{\bot}$, the canonical interval
extension ${{\boldsymbol{I}}f}$ defined in (6) is the maximal extension of $f$
among all the extensions in the continuous domain
${\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{m}_{\bot}$.
In particular, ${{\boldsymbol{I}}f}$ is Scott-continuous.

###### Proof.

Given a map $f:\mathbb{R}^{n}\to\mathbb{K}\mathbb{R}^{m}_{\bot}$, we define
$f^{\Box}:\mathbb{R}^{n}\to{\mathbb{I}}\mathbb{R}^{m}_{\bot}$ by
$f^{\Box}\coloneqq(\cdot)^{\Box}\circ f$, i. e., $\forall
x\in\mathbb{R}^{n}:f^{\Box}(x)=f(x)^{\Box}$. If $f$ is Euclidean-Scott-
continuous, then, by Corollary 2.17, so is $f^{\Box}$. It is straightforward
to verify that a map
$u:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{m}_{\bot}$
is an interval approximation of $f$ if and only if it is an interval
approximation of $f^{\Box}$. Thus, it suffices to prove the proposition for
the special case of $f:\mathbb{R}^{n}\to{\mathbb{I}}\mathbb{R}^{m}_{\bot}$.
This has been proved in [28, Lemma 3.4] for $n=m=1$. The proof given in [28],
however, is independent of the values of $n$ and $m$, and the main property
that is required is that ${\mathbb{I}}\mathbb{R}^{m}_{\bot}$ is a continuous
$\bigsqcap$-semilattice, for any $m\in\mathbb{N}$. ∎

As the restriction of the Scott topology of $\mathbb{K}\mathbb{R}^{m}_{\bot}$
over $\mathbb{R}^{m}$ is the Euclidean topology, we may consider any
continuous map $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ also as a function of type
$\mathbb{R}^{n}\to\mathbb{K}\mathbb{R}^{m}_{\bot}$, and construct its
canonical interval extension accordingly, which will be Scott-continuous.

######  Remark 2.21.

The distinction between extension and approximation is crucial in the current
article.

##  3 Global Robustness Analysis

A domain-theoretic framework for abstract robustness analysis was introduced
by [60]. The essence of the framework is in the link between robustness—a
concept defined based on perturbations—and Scott topology—a concept defined in
terms of order theory.

Assume that $\mathbb{S}$ is a metric space, $\mathbb{K}\mathbb{S}$ is the
PDCPO of non-empty compact subsets of $\mathbb{S}$, ordered under reverse
inclusion, with $\mathbb{S}$ added as the bottom element if $\mathbb{S}$ is
not compact. When $\mathbb{S}$ is a compact metric space—which is the relevant
case in our discussion—the PDCPO $\mathbb{K}\mathbb{S}$ becomes an
$\omega$-continuous domain [20, Proposition 3.4]. For each
$X\subseteq\mathbb{S}$, let
$\uparrow\negthickspace{X}\coloneqq\\{{C\in\mathbb{K}\mathbb{S}}\mid{C\subseteq
X}\\}$. The robust topology is defined over $\mathbb{K}\mathbb{S}$.
Intuitively, we regard a collection $U\subseteq\mathbb{K}\mathbb{S}$ of
compact subsets of $\mathbb{S}$ as robust open if $U$ is closed under
sufficiently small perturbations. Formally, we say that
$U\subseteq\mathbb{K}\mathbb{S}$ is _robust open_ iff $\forall C\in
U:\exists\delta>0:\uparrow\negthickspace{B(C,\delta)}\subseteq U$, in which
$B(C,\delta)$ is as defined in (2). This topology indeed captures the notion
of robustness [60, Theorem A.2]. Thus, we define:

######  Definition 3.1 (Robust Map).

Assume that $\mathbb{S}_{1}$ and $\mathbb{S}_{2}$ are compact metric spaces.
We say that $f:\mathbb{K}\mathbb{S}_{1}\to\mathbb{K}\mathbb{S}_{2}$ is robust
iff it is continuous with respect to the robust topologies on
$\mathbb{K}\mathbb{S}_{1}$ and $\mathbb{K}\mathbb{S}_{2}$.

######  Theorem 3.2 ([60, Theorem A.4]).

For any compact metric space $\mathbb{S}$, the Scott and robust topologies on
$\mathbb{K}\mathbb{S}$ coincide.

Assume that $N:\mathbb{S}_{1}\to\mathbb{S}_{2}$ is a network, and define the
reachability map $A_{N}:\mathbb{K}\mathbb{S}_{1}\to\mathbb{K}\mathbb{S}_{2}$
as follows:

$\forall X\in\mathbb{K}\mathbb{S}_{1}:\quad
A_{N}(X)\coloneqq\bigcap\\{{C\in\mathbb{K}\mathbb{S}_{2}}\mid{N(X)\subseteq
C}\\}.$ (7)

The map $A_{N}$ is clearly monotonic, that is, $\forall
X,Y\in\mathbb{K}\mathbb{S}_{1}:X\subseteq Y\implies A_{N}(X)\subseteq
A_{N}(Y)$. In case $A_{N}$ is not robust, a valid question is whether it is
possible to approximate $A_{N}$ with a robust map without losing too much
accuracy. For compact metric spaces, the answer is affirmative:

######  Corollary 3.3.

If $\mathbb{S}_{1}$ and $\mathbb{S}_{2}$ are compact metric spaces, then any
monotonic map $f:\mathbb{K}\mathbb{S}_{1}\to\mathbb{K}\mathbb{S}_{2}$ has a
‘tightest’ robust approximation
$f^{\circledcirc}:\mathbb{K}\mathbb{S}_{1}\to\mathbb{K}\mathbb{S}_{2}$
satisfying:

$\forall X\in\mathbb{K}\mathbb{S}_{1}:\quad
f^{\circledcirc}(X)=\bigsqcup\\{{f(b)}\mid{b\ll X}\\},$ (8)

which by Theorem 3.2, must be Scott-continuous. Here, by the tightest we mean
the largest under the pointwise ordering of functions.

###### Proof.

The existence of a tightest robust approximation follows from [60, Corollaries
4.4 and 4.5]. The formulation (8) follows from [60, Theorem 5.20]. ∎

In a typical application, the input space of a neural network is a bounded
region of some Euclidean space. Therefore, we consider feedforward neural
networks with input space $[-M,M]^{n}$, for some $M>0$. For the output, we
will mainly focus on the regressors of type $N:[-M,M]^{n}\to\mathbb{R}$. The
reason is that, within the framework that we have adopted, robustness analysis
of classifiers reduces to Lipschitz estimation of regressors. Assume that
$\hat{N}:[-M,M]^{n}\to\left\\{{c_{1},\ldots,c_{k}}\right\\}$ is a classifier
network. A common architecture for such a classifier comprises a regressor
$N\equiv(N_{1},\ldots,N_{k}):[-M,M]^{n}\to\mathbb{R}^{k}$ followed by an
$\arg\max$ operation. Thus, for a given input $x_{0}\in[-M,M]^{n}$, the index
of the class predicted by $\hat{N}$ is $i_{0}=\arg\max_{1\leq i\leq
k}N_{i}(x_{0})$. In certified local robustness analysis, the aim is to provide
a radius $\beta>0$ such that, for any point $x$ in the neighborhood of radius
$\beta$ around $x_{0}$, we have $\hat{N}(x)=\hat{N}(x_{0})$, i. e., a
neighborhood within which no adversarial examples can be found. A lower bound
for such a radius can be obtained by estimating the Lipschitz constants of the
functions $\\{{N_{i_{0}}-N_{i}}\mid{1\leq i\leq k,i\neq i_{0}}\\}$. As such,
robustness analysis of classifiers reduces to Lipschitz estimation of
regressors with one output. For more details, the reader may refer to [85,
Section 3].

When $N:[-M,M]^{n}\to\mathbb{R}$ is a continuous function—which is the
relevant case in our discussion—the image of the compact set $[-M,M]^{n}$ will
also be compact, which we may assume to be included in a compact interval
$[-M^{\prime},M^{\prime}]$, for $M^{\prime}\geq 0$ large enough. Furthermore,
the reachability map $A_{N}$ of (7) satisfies $\forall
X\in\mathbb{K}[-M,M]^{n}:A_{N}(X)=N(X)$.

Before focusing on regressors, however, we briefly discuss some issues related
to robustness analysis of classifiers.

###  3.1 Classifiers

Let $C=\left\\{{c_{1},\ldots,c_{k}}\right\\}$ be a set of classes, with $k\geq
1$, and assume that the function $N:[-M,M]^{n}\to C$ is (the semantic model
of) a classifier. As $C$ is discrete, it is natural to metrize $C$ using the
discrete metric:

$\forall x,y\in C:\quad d(x,y)\coloneqq\left\\{\begin{array}[]{ll}0,&\text{if
}x=y,\\\ 1,&\text{if }x\neq y.\\\ \end{array}\right.$

The main problem with a classifier of type $N:[-M,M]^{n}\to C$ is the mapping
of a connected space to a discrete one, where the topology on $[-M,M]^{n}$ is
the Euclidean topology, and the topology on $C$ is the discrete topology:

######  Proposition 3.4.

A classifier $N:[-M,M]^{n}\to C$ is continuous iff $N$ is constant.

###### Proof.

The $\Leftarrow$ direction is obvious. For the $\Rightarrow$ direction, as $C$
is a discrete space, then each singleton $\left\\{{c_{i}}\right\\}$ is an open
set. If $N$ is continuous, then each inverse image $N^{-1}(c_{i})$ must also
be open. If $N$ takes more than one value, then $[-M,M]^{n}=\bigcup_{1\leq
i\leq k}N^{-1}(c_{i})$ must be the union of at least two disjoint non-empty
open sets, which is impossible as $[-M,M]^{n}$ is a connected space. ∎

Proposition 3.4 implies that a classifier may not be computable—in the
framework of TTE [84]—unless if it is constant. This is because, by [53,
Theorem 2.13], if $N$ is computable, then it must be continuous.

By moving to compact subsets, however, it is possible to obtain tight robust
approximations of classifiers, even in the presence of discontinuities. As $C$
is a finite set, it is compact under the discrete topology, and $\mathbb{K}C$
is the set of non-empty subsets of $C$, ordered under reverse inclusion.

######  Theorem 3.5.

Assume that $A_{N}:\mathbb{K}[-M,M]^{n}\to\mathbb{K}C$ is the reachability map
for a classifier $N:[-M,M]^{n}\to C$, defined in (7). Then, $A_{N}$ has a
tightest robust approximation
$A_{N}^{\circledcirc}:\mathbb{K}[-M,M]^{n}\to\mathbb{K}C$ which is also Scott-
continuous.

###### Proof.

Since both $[-M,M]^{n}$ and $C$ are compact, the result follows from Corollary
3.3. ∎

######  Example 3.6.

Assume that the classifier $N:[-1,1]\to\left\\{{c_{1},c_{2}}\right\\}$
satisfies:

$\forall x\in[-1,1]:\quad N(x)=\left\\{\begin{array}[]{ll}c_{1},&\text{if
}x\in[-1,0],\\\ c_{2},&\text{if }x\in(0,1].\\\ \end{array}\right.$

Then, according to (7), we have:

$\forall x\in[-1,1]:\quad
A_{N}(\left\\{{x}\right\\})=\left\\{\begin{array}[]{ll}\left\\{{c_{1}}\right\\},&\text{if
}x\in[-1,0],\\\ \left\\{{c_{2}}\right\\},&\text{if }x\in(0,1],\\\
\end{array}\right.$

while according to (8), we must have:

$\forall x\in[-1,1]:\quad
A_{N}^{\circledcirc}(\left\\{{x}\right\\})=\left\\{\begin{array}[]{ll}\left\\{{c_{1}}\right\\},&\text{if
}x\in[-1,0),\\\ \left\\{{c_{1},c_{2}}\right\\},&\text{if }x=0,\\\
\left\\{{c_{2}}\right\\},&\text{if }x\in(0,1].\\\ \end{array}\right.$

As can be seen from Example 3.6, the tightest robust approximation of a non-
constant classifier may be _multi-valued_ over certain parts of the input
domain.

######  Theorem 3.7.

Assume that $N:[-M,M]^{n}\to C$ is a non-constant classifier, and let
$\hat{N}:\mathbb{K}[-M,M]^{n}\to\mathbb{K}C$ be any robust approximation of
the reachability map $A_{N}$. Then, $\hat{N}(\left\\{{x}\right\\})$ must be
multivalued for some $x\in[-M,M]^{n}$.

###### Proof.

To obtain a contradiction, assume that $\hat{N}(x)$ is single-valued for all
$x\in[-M,M]^{n}$. Therefore, the restriction of $\hat{N}$ over
$[-M,M]^{n}$—which is the set of maximal elements of
$\mathbb{K}[-M,M]^{n}$—must be equal to $N$. By Corollary 2.11, the classifier
$N:[-M,M]^{n}\to C$ must be continuous, which, combined with Proposition 3.4,
implies that $N$ is a constant classifier, which is a contradiction. ∎

In Theorem 3.7, assume that $C=\left\\{{c_{1},\ldots,c_{k}}\right\\}$. In case
$\hat{N}$ is the _tightest_ robust approximation $A_{N}^{\circledcirc}$ of
$A_{N}$, then $\hat{N}$ is required to be multivalued only over the boundaries
of $N^{-1}(c_{1}),\ldots,N^{-1}(c_{k})$. In typical applications, these
boundaries have Lebesgue measure zero, and the restriction of
$A_{N}^{\circledcirc}$ over $[-M,M]^{n}$ is single-valued and equal to $N$,
_almost everywhere_ , e. g., as in Example 3.6. Computability of
$A_{N}^{\circledcirc}$, however, is an open problem in general, and even in
cases where it is known to be computable, little is known about complexity of
its computation. Thus, in practice, one aims for robust approximations
$\hat{N}$ that are ‘sufficiently’ tight, and over $[-M,M]^{n}$, may be
multivalued on sets that are not Lebesgue measure zero.

###  3.2 Regressors

Assume that a given feedforward regressor $N:[-M,M]^{n}\to\mathbb{R}$ has $H$
hidden layers. We let $n_{0}\coloneqq n$, and for each $1\leq i\leq H$, let
$n_{i}$ denote the number of neurons in the $i$-th layer. Then, for each
$1\leq i\leq H$, we represent:

  * •

the weights as a matrix $W_{i}$, with $\dim(W_{i})=n_{i}\times n_{i-1}$;

  * •

the biases as a vector $b_{i}$, with $\dim(b_{i})=n_{i}\times 1$;

  * •

and the activation function (seen as a vector field) by
$\sigma_{i}:\mathbb{R}^{n_{i}}\to\mathbb{R}^{n_{i}}$.

From the last hidden layer to the single output neuron, the weights are
represented as the vector $W_{H+1}$, with $\dim(W_{H+1})=1\times n_{H}$, and
there are no biases or activation functions. For any input $x\in[-M,M]^{n}$,
the output of the neural network may be calculated layer-by-layer as follows:

$\left\\{\begin{array}[]{ll}Z_{0}(x)=x,\\\
\hat{Z}_{i}(x)=W_{i}Z_{i-1}(x)+b_{i},&(1\leq i\leq H),\\\
Z_{i}(x)=\sigma_{i}\left(\hat{Z}_{i}(x)\right),&(1\leq i\leq H),\\\
N(x)=W_{H+1}Z_{H}(x).\end{array}\right.$ (9)

We let ${\cal F}$ denote the class of feedforward networks with a single
output neuron, and Lipschitz continuous activation functions, whose output may
be described using (9). Each $N\in{\cal F}$ is continuous. Hence, it maps the
compact set $[-M,M]^{n}$ to a compact set, which we assume is included in
$[-M^{\prime},M^{\prime}]$, for $M^{\prime}$ large enough.

######  Proposition 3.8.

Assume that $N:[-M,M]^{n}\to[-M^{\prime},M^{\prime}]$ is a regressor in ${\cal
F}$, and $A_{N}$ is its reachability map. Then, the map $A_{N}$ has a tightest
robust approximation which is a robust extension of $N$, i. e., an extension
according to Definition 2.19 (i).

###### Proof.

Since both $[-M,M]^{n}$ and $[-M^{\prime},M^{\prime}]$ are compact, by
Corollary 3.3, the reachability map $A_{N}$ has a tightest robust
approximation
$A_{N}^{\circledcirc}:\mathbb{K}[-M,M]^{n}\to\mathbb{K}[-M^{\prime},M^{\prime}]$.
To prove that $A_{N}^{\circledcirc}$ is indeed an extension of $N$, we must
show that over the maximal elements of $\mathbb{K}[-M,M]^{n}$, the map
$A_{N}^{\circledcirc}$ takes singleton values. This is a straightforward
consequence of the fact that $N$ is continuous. ∎

This is in contrast with the case of classifiers, where the network itself did
not have to be even continuous, and we had to be content with tightest robust
approximations.

Of particular interest is the interval enclosure of a given regressor, which
may not be the tightest robust extension of the regressor, but it is a robust
extension nonetheless, which is commonly used in interval analysis of neural
networks:

######  Theorem 3.9.

Let $N:[-M,M]^{n}\to[-M^{\prime},M^{\prime}]$ be a regressor in ${\cal F}$,
and assume that:

  * •

$(\cdot)^{\Box}:\mathbb{K}[-M,M]^{n}\to{\mathbb{I}}[-M,M]^{n}$ is the hyperbox
closure map of (4).

  * •

$\tilde{N}:{\mathbb{I}}[-M,M]^{n}\to{\mathbb{I}}[-M^{\prime},M^{\prime}]$ is
any Scott-continuous extension of $N$, e. g., the canonical interval extension
of Definition 2.19 (iii).

  * •

$\iota:{\mathbb{I}}[-M^{\prime},M^{\prime}]\to\mathbb{K}[-M^{\prime},M^{\prime}]$
is the embedding of the subdomain ${\mathbb{I}}[-M^{\prime},M^{\prime}]$ into
$\mathbb{K}[-M^{\prime},M^{\prime}]$, as in Corollary 2.18.

Then, the map
$\hat{N}:\mathbb{K}[-M,M]^{n}\to\mathbb{K}[-M^{\prime},M^{\prime}]$ defined by
$\forall
X\in\mathbb{K}[-M,M]^{n}:\hat{N}(X)\coloneqq\iota(\tilde{N}(X^{\Box}))$ is a
robust extension of $N$.

###### Proof.

We know that $(\cdot)^{\Box}$, $\tilde{N}$, and $\iota$ are Scott-continuous.
Therefore, the map $\hat{N}$—being a composition of the three—is Scott-
continuous. Since both $[-M,M]^{n}$ and $[-M^{\prime},M^{\prime}]$ are
compact, by Theorem 3.2, the map $\hat{N}$ is robust. Hence, $\hat{N}$ is a
robust approximation of $N$. To prove that $\hat{N}$ is a robust extension of
$N$, we must show that over maximal input values it returns singletons. But
this is again a straightforward consequence of continuity of $N$. ∎

##  4 Validated Local Robustness Analysis

The results of Section 3 concern globally robust approximation of neural
networks. In the machine learning literature, robustness analysis is mainly
carried out locally, that is, over specific locations in the input domain. A
classifier $N:\mathbb{R}^{n}\to\left\\{{c_{1},\ldots,c_{k}}\right\\}$ may be
regarded as robust at $v\in\mathbb{R}^{n}$ if, for some neighborhood
$B(\left\\{{v}\right\\},r_{0})$ of radius $r_{0}>0$, we have: $\forall x\in
B(\left\\{{v}\right\\},r_{0}):N(x)=N(v)$. In the current article, by attack-
agnostic robustness measurement we mean obtaining a tight lower bound for the
largest such $r_{0}$. This, in turn, can be reduced to obtaining a tight upper
bound for the local Lipschitz constant in a neighborhood of the point $v$ [45,
85], which is the focus of the rest of the current article.

From a theoretical angle, the more challenging networks for Lipschitz analysis
are those with non-differentiable activation functions, such as
$\operatorname{\mathrm{ReLU}}$. For such functions, a generalized notion of
gradient must be used, e. g., Clarke-gradient [14]. Indeed, Clarke-gradient
has been used recently for robustness analysis of Lipschitz neural networks
[49, 50]. On the other hand, the _domain-theoretic $L$-derivative_ (Definition
4.4 LABEL:[below]def:L_derivative) introduced by [25] also coincides with
Clarke-gradient, a fact which was proven for finite-dimensional Banach spaces
by [25], and later generalized to infinite-dimensional Banach spaces by [46].
Using this, we develop a domain-theoretic framework for Lipschitz analysis of
feedforward networks, which provides a theoretical foundation for methods such
as [49, 50].

Working with the $L$-derivative, however, is computationally costly. It
requires computation over $\mathbb{C}\mathbb{R}^{n}_{\bot}$ and involves a
rather complicated chain rule (Proposition 4.3). Therefore, we consider a
hyperbox approximation of the $L$-derivative, which we call the
$\hat{L}$-derivative (Definition 4.7). This will allow us to work with the
relatively simpler domain ${\mathbb{I}}\mathbb{R}^{n}_{\bot}$ and also a much
simpler chain rule (Lemma 4.9). Finally, we will move on to computation of the
Lipschitz constant of a given regressor, according to the following overall
strategy.

Assume that $N\in{\cal F}$ is a feedforward regressor. By using (9), we obtain
a unique closed-form expression for $N$ in terms of the weights, biases, and
the activation functions. Let
$\tilde{L}(N):{\mathbb{I}}[-M,M]^{n}\to{\mathbb{I}}\mathbb{R}^{n}_{\bot}$
denote the interval approximation of the $\hat{L}$-derivative of $N$ obtained
by applying the chain rule of Lemma 4.9 over this closed-form expression. The
map $\tilde{L}(N)$ is Scott-continuous. Since we consider
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{\infty}$ as the
perturbation norm, we must compute the
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{1}$-norm of the gradients
(see Remark 2.1). First, we extend the absolute value function to real
intervals by defining
$\lvert\thinspace{I}\thinspace\rvert\coloneqq\\{{\lvert\thinspace{x}\thinspace\rvert}\mid{x\in
I}\\}$ for any interval $I$. Then, we extend the
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{1}$ norm to hyperboxes by
defining
${\left\lVert\thinspace{\prod_{i=1}^{n}I_{i}}\thinspace\right\rVert}_{1}\coloneqq\sum_{i=1}^{n}\lvert\thinspace{I_{i}}\thinspace\rvert$.
For instance, we have
${\left\lVert\thinspace{[-1,1]\times[-1,2]}\thinspace\right\rVert}_{1}=[0,1]+[0,2]=[0,3]$.
The function
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{1}:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
is also Scott-continuous. Therefore, the composition
${\left\lVert\thinspace{\tilde{L}(N)}\thinspace\right\rVert}_{1}:{\mathbb{I}}[-M,M]^{n}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
is a Scott-continuous function. All that remains is to find the maximum value
and the maximum set of
${\left\lVert\thinspace{\tilde{L}(N)}\thinspace\right\rVert}_{1}$ to
approximate the Lipschitz constant and the set of points where the maximum
Lipschitz values are attained, within any given neighborhood. This is one of
the main contributions of the current article and will be presented in detail
in Section 4.4.

If, instead of ${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, we use the finer domain
$\mathbb{C}\mathbb{R}^{n}_{\bot}$, then the suitable chain rule will be that
of Proposition 4.3, which provides a more accurate result, at the cost of
calculating the closed convex hull. This has also been observed in [50]. By
using the $\hat{L}$-derivative and the chain rule of Lemma 4.9, we avoid
computation of the convex hull, and trade some accuracy for efficiency. As it
turns out, however, _we will not lose any accuracy over general position
$\operatorname{\mathrm{ReLU}}$ networks (Theorem 4.19), or over differentiable
networks (Theorem 4.21)_.

######  Remark 4.1.

In this article, we focus on the
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{\infty}$-norm for
perturbations because we will be working with hyperboxes, and in the
implementation, we will use interval arithmetic. This does not mean, however,
that we cannot analyze robustness with respect to
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{p}$ for other values of
$p$. For example, assume that we are given a closed neighborhood
$V(x_{0},r)\coloneqq\\{{x\in\mathbb{R}^{n}}\mid{{\left\lVert\thinspace{x-x_{0}}\thinspace\right\rVert}_{p}\leq
r}\\}$, for some $p\in[1,\infty)$. We can first cover the neighborhood
$V(x_{0},r)$ with (a large number of sufficiently small) hyperboxes
$B=\\{{b_{i}}\mid{i\in I}\\}$. We obtain an overapproximation of the
${\left\lVert\thinspace{.}\thinspace\right\rVert}_{p^{\prime}}$-norm of the
Clarke-gradient over each box $b_{i}$, and then obtain an approximation of the
Lipschitz constant over the neighborhood $V(x_{0},r)$.

###  4.1 Lipschitz Continuity and Lipschitz Constant

Assume that $(X,\delta_{X})$ and $(Y,\delta_{Y})$ are two metric spaces,
$f:X\to Y$, and $A\subseteq X$. A real number $L\geq 0$ is said to be a
_Lipschitz bound_ for $f$ over $A$ if:

$\forall a,b\in A:\delta_{Y}(f(a),f(b))\leq L\delta_{X}(a,b).$ (10)

A function for which such a bound exists is said to be Lipschitz continuous
over $A$, and the smallest $L$ satisfying (10)—which must exist—is called the
_Lipschitz constant_ of $f$ over $A$.

If $f$ is Lipschitz continuous over $A$, then it is also continuous over $A$,
although the converse is not true. For instance, the function
$h:\ell^{1}_{\infty}\to\ell^{1}_{\infty}$, defined by
$h(x)\coloneqq\sqrt{\lvert\thinspace{x}\thinspace\rvert}$, is continuous
everywhere, but not Lipschitz continuous over any neighborhood of $0$.
Whenever $X$ is a compact subset of $\ell^{n}_{\infty}$ and $f:X\to\mathbb{R}$
is continuously differentiable, then $f$ is Lipschitz over $X$, but again, the
converse is not true. For example, the function $h:[-1,1]\to\mathbb{R}$,
defined by $h(x)\coloneqq\lvert\thinspace{x}\thinspace\rvert$, is Lipschitz
continuous, but not differentiable at $0$. To summarize, we have the following
inclusions for functions over compact subsets of $\mathbb{R}^{n}$:

$\text{Continuously Differentiable}\subset\text{Lipschitz
Continuous}\subset\text{Continuous}.$

Even though the inclusions are strict, there is a close relationship between
Lipschitz continuity and differentiability. Rademacher’s theorem states that
if $f$ is Lipschitz continuous over an open subset $A$ of $\mathbb{R}^{n}$,
then it is (Fréchet) differentiable almost everywhere (with respect to the
Lebesgue measure) over $A$ [13, Corollary 4.19]. For example, the function
$\operatorname{\mathrm{ReLU}}$—defined by $\forall
x\in\mathbb{R}:\operatorname{\mathrm{ReLU}}(x)\coloneqq\max(0,x)$—is
(classically) differentiable everywhere except at zero:

$(\operatorname{\mathrm{ReLU}})^{\prime}(x)=\left\\{\begin{array}[]{ll}1,&\text{if
}x>0,\\\ 0,&\text{if }x<0,\\\ \mathrm{undefined},&\text{if
}x=0.\end{array}\right.$

###  4.2 Clarke-Gradient

Although $\operatorname{\mathrm{ReLU}}$ is not differentiable at $0$, one-
sided limits of the derivative exist, i. e.:

$\lim_{x\to
0^{-}}(\operatorname{\mathrm{ReLU}})^{\prime}(x)=0\quad\text{and}\quad\lim_{x\to
0^{+}}(\operatorname{\mathrm{ReLU}})^{\prime}(x)=1.$ (11)

As shown in Figure 1 (Left), in geometric terms, for every $\lambda\in[0,1]$,
the line $y_{\lambda}(x)\coloneqq\lambda x$ satisfies:

$y_{\lambda}(0)=\operatorname{\mathrm{ReLU}}(0)\quad\wedge\quad\forall
x\in\mathbb{R}:y_{\lambda}(x)\leq\operatorname{\mathrm{ReLU}}(x).$ (12)

The interval $[0,1]$ is the largest set of $\lambda$ values which satisfy
(12), and may be considered as a gradient _set_ for
$\operatorname{\mathrm{ReLU}}$ at zero. The concept of generalized (Clarke)
gradient formalizes this idea. Clarke-gradient may be defined for functions
over any Banach space $X$ [14, page 27]. When $X$ is finite-dimensional,
however, the Clarke-gradient of $f:X\to\mathbb{R}$ has a simpler
characterization. Assume that $X$ is an open subset of $\mathbb{R}^{n}$, and
$f:X\to\mathbb{R}$ is Lipschitz continuous over $X$. By Rademacher’s theorem,
we know that the set ${\cal N}_{f}\coloneqq\\{{z\in X}\mid{f\text{ is not
differentiable at }z\,}\\}$ has Lebesgue measure zero. The Clarke-gradient
$\partial f(x)$ of $f$ at any $x\in X$ satisfies:

$\partial
f(x)=\operatorname{co}\left\\{{\lim_{i\to\infty}f^{\prime}(x_{i})}\;\vrule\;{\lim_{i\to\infty}x_{i}=x\wedge\forall
i\in\mathbb{N}:x_{i}\notin{\cal N}_{f}}\right\\}.$ (13)

Equation (13) should be interpreted as follows: take any sequence
$(x_{i})_{i\in\mathbb{N}}$, which completely avoids the set ${\cal N}_{f}$ of
points where $f$ is not differentiable, converges to $x$, and for which
$\lim_{i\to\infty}f^{\prime}(x_{i})$ exists. The Clarke-gradient is the convex
hull of these limits [14, page 63]. It should be clear from (11) and (13) that
$\partial\operatorname{\mathrm{ReLU}}(0)=[0,1]$.

######  Proposition 4.2 ([14, Proposition 2.1.2]).

Assume that $p\in[1,\infty]$. Let $X$ be an open subset of $\ell^{n}_{p}$ and
$f:X\to\mathbb{R}$ be Lipschitz continuous over $X$ with a Lipschitz bound
$L$. Then, for every $x\in X$:

  1. (i)

$\partial f(x)$ is a non-empty, convex, and compact subset of
$\mathbb{R}^{n}$.

  2. (ii)

$\forall\xi\in\partial
f(x):{\left\lVert\thinspace{\xi}\thinspace\right\rVert}_{p^{\prime}}\leq L$,
in which $p^{\prime}$ is the conjugate of $p$ as defined in (1).

The following chain rule for the Clarke-gradient is crucial for proving the
soundness of methods such as [11, 10, 49, 50, 57]:

######  Proposition 4.3 (Chain Rule, [14, Theorem 2.3.9]).

Assume that $X\subseteq\mathbb{R}^{n}$ is open,
$h=(h_{1},\ldots,h_{m}):X\to\mathbb{R}^{m}$ is Lipschitz continuous in a
neighborhood of $x\in X$, and $g:\mathbb{R}^{m}\to\mathbb{R}$ is Lipschitz
continuous in a neighborhood of $h(x)\in\mathbb{R}^{m}$. Then, $f\coloneqq
g\circ h$ is Lipschitz in a neighborhood of $x$, and we have:

$\partial
f(x)\subseteq\overline{\operatorname{co}}\left\\{{\sum_{i=1}^{m}\alpha_{i}\zeta_{i}}\;\vrule\;{\zeta_{i}\in\partial
h_{i}(x)\wedge\alpha\in\partial g(h(x))}\right\\},$ (14)

in which, $\overline{\operatorname{co}}$ denotes the closed convex hull.

Unlike the chain rule for the classical derivative, the chain rule of (14)
provides only a subset relation, which, in general, is not an equality. A
simple example of lack of equality is shown in Figure 2, which demonstrates an
instance of the dependency problem (see Section 4.6.2).

###  4.3 Domain-Theoretic $\boldsymbol{L}$-Derivative

We recall the concept of a domain-theoretic derivative for a function
$f:U\to\mathbb{R}$ defined on an open set $U\subseteq\mathbb{R}^{n}$. Assume
that $(X,\Omega)$ is a topological space, and $(D,\sqsubseteq)$ is a PDCPO,
with bottom element $\bot$. Then, for any open set $O\in\Omega$, and any
element $b\in D$, we define the single-step function $b\chi_{O}:X\to D$ as
follows:

$b\chi_{O}(x)\coloneqq\left\\{\begin{array}[]{ll}b,&\text{if }x\in O,\\\
\bot,&\text{if }x\in X\setminus O.\end{array}\right.$ (15)

Figure 1: Left: The function
$\operatorname{\mathrm{ReLU}}:(-2,2)\to\mathbb{R}$ (solid black) and the line
$y_{\lambda}$ (in blue) from (12) for some $\lambda\in[0,1]$, which lies
entirely in the shaded cone. Right: The step function
$b\chi_{O}:(-2,2)\to{\mathbb{I}}\mathbb{R}_{\bot}$ in which $b=[0,1]$ and
$O=(-1,1)$. We have $\operatorname{\mathrm{ReLU}}\in\delta(O,b)$.

######  Definition 4.4 ($L$-derivative [21, Section 2]).

Assume that $O\subseteq U\subseteq\mathbb{R}^{n}$, both $O$ and $U$ are open,
and $b\in\mathbb{C}\mathbb{R}^{n}_{\bot}$:

  1. (i)

The single-step tie $\delta(O,b)$ is the set of all functions
$f:U\to\mathbb{R}$ that satisfy $\forall x,y\in O:b(x-y)\sqsubseteq
f(x)-f(y)$. The set $b(x-y)$ is obtained by taking the inner product of every
element of $b$ with the vector $x-y$, and $\sqsubseteq$ is the reverse
inclusion order on $\mathbb{C}\mathbb{R}^{1}_{\bot}$.

  2. (ii)

The $L$-derivative of any function $f:U\to\mathbb{R}$ is defined as:

$L(f)\coloneqq\bigsqcup\left\\{{b\chi_{O}}\;\vrule\;{O\subseteq
U,b\in\mathbb{C}\mathbb{R}^{n}_{\bot},f\in\delta(O,b)}\right\\}.$

As such, the map $L$ has type
$L:(U\to\mathbb{R})\to(U\to\mathbb{C}\mathbb{R}^{n}_{\bot})$.

######  Example 4.5.

Assume that $U=(-2,2)$. Consider the function
$\operatorname{\mathrm{ReLU}}:(-2,2)\to\mathbb{R}$ and the step function
$b\chi_{O}:(-2,2)\to{\mathbb{I}}\mathbb{R}_{\bot}$ with $b=[0,1]$ and
$O=(-1,1)$, as depicted in Figure 1. We have
$\operatorname{\mathrm{ReLU}}\in\delta(O,b)$. In fact, we have
$\operatorname{\mathrm{ReLU}}\in\delta(O,b)$ for any open set $O$ containing
zero, e. g., the interval $(-2^{n},2^{n})$, for any $n\in\mathbb{N}$. Thus,
one may verify that:

$\forall x\in(-2,2):\quad
L(\operatorname{\mathrm{ReLU}})(x)=\left\\{\begin{array}[]{ll}0,&\text{if
}x<0,\\\ {[0,1]},&\text{if }x=0,\\\ 1,&\text{if }x>0.\end{array}\right.$

For any $f:U\to\mathbb{R}$, the $L$-derivative
$L(f):U\to\mathbb{C}\mathbb{R}^{n}_{\bot}$ is Euclidean-Scott-continuous. When
$f$ is classically differentiable at $x\in U$, the $L$-derivative and the
classical derivative coincide [25]. Many of the fundamental properties of the
classical derivative can be generalized to the domain-theoretic one, e. g.,
additivity and the chain rule [23]. Of particular importance to the current
paper is the following result:

######  Theorem 4.6 ([25, 46]).

Over Lipschitz continuous functions $f:U\to\mathbb{R}$, the domain-theoretic
$L$-derivative coincides with the Clarke-gradient, i. e., $\forall x\in
U:L(f)(x)=\partial f(x)$.

In this paper, instead of working with the general convex sets in
$\mathbb{C}\mathbb{R}^{n}_{\bot}$, we work with the simpler hyperboxes in
${\mathbb{I}}\mathbb{R}^{n}_{\bot}$:

######  Definition 4.7 ($\hat{L}$-derivative).

Let $U\subseteq\mathbb{R}^{n}$ be an open set. In the definition of
$L$-derivative (Definition 4.4), if $\mathbb{C}\mathbb{R}^{n}_{\bot}$ is
replaced with ${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, then we obtain the concept
of $\hat{L}$-derivative for functions of type $U\to\mathbb{R}$. For vector
valued functions $f=(f_{1},\ldots,f_{m}):U\to\mathbb{R}^{m}$, we define:

$\hat{L}(f)\coloneqq\left(\hat{L}(f_{1}),\ldots,\hat{L}(f_{m})\right)^{\intercal},$

in which, $(\cdot)^{\intercal}$ denotes the transpose of a matrix. In other
words, for each $1\leq i\leq m$ and $x\in U$, let
$\hat{L}(f_{i})(x)\equiv(\alpha_{i,1},\ldots,\alpha_{i,n})\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}$.
Then, for each $x\in U$, $\hat{L}(f)(x)$ is the $m\times n$ interval matrix
$[\alpha_{i,j}]_{1\leq i\leq m,1\leq j\leq n}$.

The $\hat{L}$-derivative is, in general, coarser than the $L$-derivative, as
the following example demonstrates:

######  Example 4.8.

Let $D_{n}$ denote the closed unit disc in $\mathbb{R}^{n}$, and define
$f:\mathbb{R}^{n}\to\mathbb{R}$ by
$f(x)={\left\lVert\thinspace{x}\thinspace\right\rVert}_{2}$, in which
${\left\lVert\thinspace{\cdot}\thinspace\right\rVert}_{2}$ is the Euclidean
norm. Then, $L(f)(0)=D_{n}$, while $\hat{L}(f)(0)=[-1,1]^{n}$.

As such, we lose some precision by using hyperboxes. On the other hand, on
hyperboxes, the chain rule has a simpler form:

######  Lemma 4.9 (Chain Rule on Hyperboxes, [23, Lemma 3.3]).

For any two functions $g:\mathbb{R}^{k}\to\mathbb{R}^{m}$ and
$f:\mathbb{R}^{m}\to\mathbb{R}^{n}$, we have:

$\hat{L}(f)(g(x))\times\hat{L}(g)(x)\sqsubseteq\hat{L}(f\circ g)(x).$

###  4.4 Maximization Algorithm

For any closed compact interval $I$, we write
$I=[{\underline{I}},{\overline{I}}]$. Assume that
$f:\mathbb{R}^{n}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$ and
$X\subseteq\mathbb{R}^{n}$ is a compact set. We define the maximum value and
maximum set of $f$ over $X$ as follows:

$\left\\{\begin{array}[]{l}\operatorname{max_{V}}(f,X)\coloneqq[{\underline{z}},{\overline{z}}],\quad\text{where
${\underline{z}}\coloneqq\sup_{x\in X}{\underline{f(x)}}$ and
${\overline{z}}\coloneqq\sup_{x\in X}{\overline{f(x)}}$ },\\\
\operatorname{max_{S}}(f,X)\coloneqq\left\\{{x\in
X}\;\vrule\;{f(x)\cap\operatorname{max_{V}}(f,X)\neq\emptyset}\right\\}.\end{array}\right.$
(16)

For any
$\hat{f}:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$,
the maximum value and maximum set are defined by considering the restriction
of $\hat{f}$ over the maximal elements $\mathbb{R}^{n}$.

######  Proposition 4.10.

Assume that
$\hat{f}:{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
is Scott-continuous, and $X$ is a compact subset of $\mathbb{R}^{n}$. Then
$\operatorname{max_{S}}(\hat{f},X)$ is a compact subset of $\mathbb{R}^{n}$.

###### Proof.

Let $f:\mathbb{R}^{n}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$ be the restriction
of $\hat{f}$ over the maximal elements $\mathbb{R}^{n}$. As $X$ is compact, it
suffices to prove that $\operatorname{max_{S}}(f,X)$ is a closed subset of
$X$. Let $x\in X\setminus\operatorname{max_{S}}(f,X)$ be an arbitrary point,
and write $f(x)=[{\underline{f(x)}},{\overline{f(x)}}]$. Since
$x\not\in\operatorname{max_{S}}(f,X)$, then
${\overline{f(x)}}<{\underline{z}}$, where
$z\equiv[{\underline{z}},{\overline{z}}]=\operatorname{max_{V}}(f,X)$. As
$\hat{f}$ is Scott-continuous, then, by Corollary 2.11, the map $f$ must be
Euclidean-Scott-continuous. In turn, by Proposition 2.14, the function
${\overline{f}}$ must be upper semicontinuous. As
${\overline{f(x)}}<{\underline{z}}$, for some neighborhood $O_{x}$ of $x$, we
must have: $\forall y\in O_{x}:{\overline{f(y)}}<{\underline{z}}$. This
entails that $O_{x}\cap\operatorname{max_{S}}(f,X)=\emptyset$. ∎

Thus, in the sequel, we consider the functionals:

$\left\\{\begin{array}[]{ll}\operatorname{max_{V}}:&({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to\mathbb{K}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot},\\\
\operatorname{max_{S}}:&({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to\mathbb{K}\mathbb{R}^{n}_{\bot}\to\mathbb{K}\mathbb{R}^{n}_{\bot},\end{array}\right.$

which, given any Scott-continuous
$f\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
and compact set $X\in\mathbb{K}\mathbb{R}^{n}_{\bot}$, return the maximum
value and the maximum set of $f$ over $X$, respectively. As an example, assume
that
$f\in{\mathbb{I}}\mathbb{R}^{1}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
is the canonical interval extension of
$\partial\operatorname{\mathrm{ReLU}}:\mathbb{R}\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$.
Then:

$\left\\{\begin{array}[]{ll}\operatorname{max_{V}}(f,[-1,1])=[1,1],&\operatorname{max_{S}}(f,[-1,1])=[0,1],\\\
\operatorname{max_{V}}(f,[-1,0])=[0,1],&\operatorname{max_{S}}(f,[-1,0])=[-1,0],\\\
\end{array}\right.$ (17)

which shows that the functionals $\operatorname{max_{V}}$ and
$\operatorname{max_{S}}$ are not even monotone in their second argument. For
fixed second arguments, however, both are Scott-continuous in their first
argument:

######  Theorem 4.11 (Scott-continuity).

For any $X\in\mathbb{K}\mathbb{R}^{n}_{\bot}$, the functionals
$\operatorname{max_{V}}(\cdot,X):({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
and
$\operatorname{max_{S}}(\cdot,X):({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to\mathbb{K}\mathbb{R}^{n}_{\bot}$
are Scott-continuous.

###### Proof.

As the domains
${\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$,
${\mathbb{I}}\mathbb{R}^{1}_{\bot}$, and $\mathbb{K}\mathbb{R}^{n}_{\bot}$ are
$\omega$-continuous, by Proposition 2.9, it suffices to show that both
functionals are monotone and preserve the suprema of $\omega$-chains.

To prove monotonicity, assume that
$f,g\in{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$,
and $f\sqsubseteq g$. Thus, $\forall x\in
X:({\underline{f(x)}}\leq{\underline{g(x)}})\wedge({\overline{g(x)}}\leq{\overline{f(x)}})$.
This implies that $\sup_{x\in X}{\underline{f(x)}}\leq\sup_{x\in
X}{\underline{g(x)}}$ and $\sup_{x\in X}{\overline{g(x)}}\leq\sup_{x\in
X}{\overline{f(x)}}$. Therefore:

$\operatorname{max_{V}}(f,X)\sqsubseteq\operatorname{max_{V}}(g,X).$ (18)

Take an arbitrary $x\in\operatorname{max_{S}}(g,X)$. By definition, we must
have $g(x)\cap\operatorname{max_{V}}(g,X)\neq\emptyset$, which, combined with
the assumption $f\sqsubseteq g$, implies
$f(x)\cap\operatorname{max_{V}}(g,X)\neq\emptyset$. This, combined with (18),
entails that $f(x)\cap\operatorname{max_{V}}(f,X)\neq\emptyset$. Hence,
$x\in\operatorname{max_{S}}(f,X)$, and we have
$\operatorname{max_{S}}(f,X)\sqsubseteq\operatorname{max_{S}}(g,X)$. The proof
of monotonicity is complete.

Next, let $(f_{n})_{n\in\mathbb{N}}$ be an $\omega$-chain in
${\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$,
with $f=\bigsqcup_{n\in\mathbb{N}}f_{n}$. As both $\operatorname{max_{V}}$ and
$\operatorname{max_{S}}$ are monotone, we have
$\bigsqcup_{n\in\mathbb{N}}\operatorname{max_{V}}(f_{n},X)\sqsubseteq\operatorname{max_{V}}(f,X)$
and
$\bigsqcup_{n\in\mathbb{N}}\operatorname{max_{S}}(f_{n},X)\sqsubseteq\operatorname{max_{S}}(f,X)$.
Thus, it only remains to prove the $\sqsupseteq$ direction in both cases. In
what follows, we let
$z\equiv[{\underline{z}},{\overline{z}}]\coloneqq\operatorname{max_{V}}(f,X)$,
and for each $n\in\mathbb{N}$,
$z_{n}\equiv[{\underline{z_{n}}},{\overline{z_{n}}}]\coloneqq\operatorname{max_{V}}(f_{n},X)$.

$\operatorname{max_{V}}$:

First, consider an arbitrary $t<{\underline{z}}$ and let
$\epsilon\coloneqq({\underline{z}}-t)/4$. By definition, $\exists x_{0}\in
X:{\underline{f(x_{0})}}>{\underline{z}}-\epsilon$. Since
$f=\bigsqcup_{n\in\mathbb{N}}f_{n}$, we must have $\exists
N_{0}\in\mathbb{N}:\forall n\geq
N_{0}:{\underline{f_{n}(x_{0})}}>{\underline{f(x_{0})}}-\epsilon>{\underline{z}}-2\epsilon>t$.
Therefore:

$\forall t\in\mathbb{R}:\quad t<{\underline{z}}\implies\exists
n\in\mathbb{N}:t<{\underline{z_{n}}}.$ (19)

Next, we consider an arbitrary $t>{\overline{z}}$ and let
$\epsilon\coloneqq(t-{\overline{z}})/4$. Since
$f=\bigsqcup_{n\in\mathbb{N}}f_{n}$, we must have:

$\forall x\in X:\exists N_{x}\in\mathbb{N}:\forall n\geq
N_{x}:\quad{\overline{f_{n}(x)}}<{\overline{f(x)}}+\epsilon.$

For each fixed $n\in\mathbb{N}$, the map ${\overline{f_{n}}}$ is upper
semicontinuous. Hence, for any $x\in X$, there exists a neighborhood $O_{x}$
of $x$ for which we have

$\forall y\in
O_{x}:\quad{\overline{f_{N_{x}}(y)}}<{\overline{f_{N_{x}}(x)}}+\epsilon<{\overline{f(x)}}+2\epsilon\leq{\overline{z}}+2\epsilon.$

As the sequence $(f_{n})_{n\in\mathbb{N}}$ is assumed to be an $\omega$-chain,
the upper bounds $({\overline{f_{n}}})_{n\in\mathbb{N}}$ are non-increasing,
and we obtain:

$\forall x\in X:\forall y\in O_{x}:\forall n\geq
N_{x}:\quad{\overline{f_{n}(y)}}\leq{\overline{z}}+2\epsilon.$

As $X$ is compact and $X\subseteq\bigcup_{x\in X}O_{x}$, then for finitely
many neighborhoods $\left\\{{O_{x_{1}},\ldots,O_{x_{m}}}\right\\}$ we must
have $X\subseteq\bigcup_{i=1}^{m}O_{x_{i}}$. Thus, if we let
$N_{0}\coloneqq\max\left\\{{N_{x_{1}},\ldots,N_{x_{m}}}\right\\}$, we must
have:

$\forall n\geq N_{0}:\forall y\in
X:\quad{\overline{f_{n}(y)}}\leq{\overline{z}}+2\epsilon,$

which entails that:

$\forall t\in\mathbb{R}:\quad{\overline{z}}<t\implies\exists
N_{0}\in\mathbb{N}:\forall n\geq
N_{0}:{\overline{z_{n}}}\leq{\overline{z}}+2\epsilon<t.$ (20)

From (19) and (20), we deduce that
$\bigsqcup_{n\in\mathbb{N}}\operatorname{max_{V}}(f_{n},X)\sqsupseteq\operatorname{max_{V}}(f,X)$.

$\operatorname{max_{S}}$:

Assume that $x\in X\setminus\operatorname{max_{S}}(f,X)$, which implies that
${\overline{f(x)}}<{\underline{z}}$. Take
$\epsilon=\left({\underline{z}}-{\overline{f(x)}}\right)/4$. Since
${\underline{z}}=\sup_{y\in X}{\underline{f(y)}}$, there exists $x_{1}\in X$
such that ${\underline{f(x_{1})}}>{\underline{z}}-\epsilon$. As
$f=\bigsqcup_{n\in\mathbb{N}}f_{n}$, there exists $N_{1}\in\mathbb{N}$, such
that:

$\forall n\geq
N_{1}:{\underline{f_{n}(x_{1})}}>{\underline{f(x_{1})}}-\epsilon>{\underline{z}}-2\epsilon.$
(21)

Furthermore, there exists $N_{2}\in\mathbb{N}$, such that:

$\forall n\geq
N_{2}:{\overline{f_{n}(x)}}<{\overline{f(x)}}+\epsilon={\underline{z}}-3\epsilon.$
(22)

Taking $N_{3}=\max\left\\{{N_{1},N_{2}}\right\\}$, from (21) and (22), we
obtain $\forall n\geq
N_{3}:{\overline{f_{n}(x)}}<{\underline{f_{n}(x_{1})}}\leq{\underline{z_{n}}}$,
which implies that $x\not\in\operatorname{max_{S}}(f_{n},X)$. Therefore,
$\bigsqcup_{n\in\mathbb{N}}\operatorname{max_{S}}(f_{n},X)\sqsupseteq\operatorname{max_{S}}(f,X)$.

The proof is complete. ∎

Algorithm 1 presents an implementation of $\operatorname{max_{V}}$ and
$\operatorname{max_{S}}$, which we use in our experiments. As pointed out
before, the set $B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$, consisting of finite
unions of rational hyperboxes, forms a countable basis for the
$\omega$-continuous domain $\mathbb{K}\mathbb{R}^{n}_{\bot}$. We have used
this fact in Algorithm 1, which outputs enclosures of the maximum set using
elements of $B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$.

Apart from the function $f$ and the neighborhood $X$, Algorithm 1 requires
three other input parameters: $m_{\iota}$, $m_{b}$, and $\delta$. These are
crucial for making sure that the algorithm stops below a reasonable threshold.
In line 9, we discard all the boxes in $B_{i}$ which cannot possibly contain
any maximum points. How many boxes are discarded depends on a number of
factors, most importantly, whether the function $f$ has a flat section around
the maximum set or not. The condition
$2^{n}\lvert\thinspace{\hat{B}_{i}}\thinspace\rvert\geq m_{b}$ in line 13 is
necessary, because in cases where flat sections exist, the number of elements
of $\hat{B}_{i}$ grows exponentially with each iteration. In the if statement
of line 10, we make sure that bisecting $\hat{B}_{i}$ will not create too many
boxes before executing the bisection command in the first place.

Algorithm 1 Enclosing $\operatorname{max_{V}}(f,X)$ and
$\operatorname{max_{S}}(f,X)$.

0: $X\subseteq\mathbb{R}^{n}$ : A finite union of hyperboxes;
$f\in{\mathbb{I}}X\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$: A Scott-
continuous interval function; $m_{\iota}\in\mathbb{N}$: maximum iterations;
$m_{b}\in\mathbb{N}$: maximum number of boxes; $\delta\geq 0$: a dyadic number
indicating target width.

0: Enclosures of $\operatorname{max_{V}}(f,X)$ and
$\operatorname{max_{S}}(f,X)$

1: $i\leftarrow 0$ // number of iterations

2: $B_{1}\leftarrow X$ // initialize the set of hyperboxes

3: repeat

4: $i\leftarrow i+1$

5: $L_{i}\coloneqq\\{{f(b)}\mid{b\in B_{i}}\\}$

6: ${\overline{z_{i}}}\leftarrow\max\\{{{\overline{I}}}\mid{I\in L_{i}}\\}$ //
the maximum upper end point

7: ${\underline{z_{i}}}\leftarrow\max\\{{{\underline{I}}}\mid{I\in L_{i}}\\}$
// the maximum lower end point

8: $I_{i}\coloneqq[{{\underline{z_{i}}},{\overline{z_{i}}}}]$ // maximum value

9: $\hat{B}_{i}\leftarrow\\{{b\in B_{i}}\mid{{\underline{z_{i}}}\in f(b)}\\}$
// retain only those boxes that could contain maximum points

10: if $(2^{n}\lvert\thinspace{\hat{B}_{i}}\thinspace\rvert<m_{b})$ then

11: $B_{i+1}\leftarrow\mathtt{bisect}(\hat{B}_{i})$ // bisect all boxes in
$\hat{B}_{i}$, along all dimensions

12: end if

13: until $(i\geq m_{\iota})$ or
$(2^{n}\lvert\thinspace{\hat{B}_{i}}\thinspace\rvert\geq m_{b})$ or
$({\overline{z_{i}}}-{\underline{z_{i}}}\leq\delta)$

14: return $(I_{i},\hat{B}_{i})$

Algorithm 1 closely follows the definition of $\operatorname{max_{V}}$ and
$\operatorname{max_{S}}$ given in (16). Hence, the proof of soundness is
straightforward:

######  Lemma 4.12 (Soundness).

For any given $f$, $X$, $m_{\iota}$, $m_{b}$, and $\delta$, Algorithm 1 halts
in finite time and returns enclosures of $\operatorname{max_{V}}(f,X)$ and
$\operatorname{max_{S}}(f,X)$.

###### Proof.

We prove the lemma by induction over $i$, that is, the number of iterations of
the main repeat-until loop. When $i=1$, we have $B_{i}=X$, which trivially
contains $\operatorname{max_{S}}(f,X)$.

For the inductive step, assume that at iteration $i$, before the start of the
loop, the finite set of hyperboxes $B_{i}$ satisfies
$\operatorname{max_{S}}(f,X)\subseteq\bigcup B_{i}$. Assume that $\tilde{f}\in
X\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$ is the restriction of the
interval function $f$ over the maximal elements $X$, which entails that $f$ is
an interval approximation of $\tilde{f}$. As a result, $I_{i}$ is an
overapproximation of $\operatorname{max_{V}}(\tilde{f},X)$, which entails that
the set $\hat{B}_{i}$—obtained in line 9—is an overapproximation of
$\operatorname{max_{S}}(\tilde{f},X)$.

The fact that the algorithm always halts in finite time follows from the
stopping condition $i\geq m_{\iota}$ of line 13. ∎

The algorithm combines the definition of $\operatorname{max_{V}}$ and
$\operatorname{max_{S}}$ given in (16) with a bisection operation that
increases the accuracy at each iteration. We prove that the result can be
obtained to any degree of accuracy. To that end, we first define the following
sequence:

$\left\\{\begin{array}[]{ll}X_{0}\coloneqq X,&\\\
X_{i+1}\coloneqq\mathtt{bisect}(X_{i}),&(\forall
i\in\mathbb{N}).\end{array}\right.$ (23)

Note that by bisecting a box $b\subseteq\mathbb{R}^{n}$ we mean bisecting
along all the $n$ coordinates, which generates $2^{n}$ boxes. It is
straightforward to verify that $\forall i\in\mathbb{N}:\hat{B}_{i}\subseteq
B_{i}\subseteq X_{i}$.

######  Theorem 4.13 (Completeness).

Assume that $\Psi(\hat{f},X,m_{\iota},m_{b},\delta)$ is the function
implemented by Algorithm 1. For any given $\hat{f}$ and $X$:

$\lim_{m_{\iota}\to\infty,m_{b}\to\infty,\delta\to
0}\Psi(\hat{f},X,m_{\iota},m_{b},\delta)=\left(\operatorname{max_{V}}(\hat{f},X),\operatorname{max_{S}}(\hat{f},X)\right).$

###### Proof.

We write the restriction of
$\hat{f}\in{\mathbb{I}}X\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$ over the
maximal elements as
$f\equiv[{\underline{f}},{\overline{f}}]:X\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$.
As $\hat{f}$ is Scott-continuous, by Corollary 2.11 and Proposition 2.14,
${\underline{f}}$ and ${\overline{f}}$ are lower and upper semicontinuous,
respectively.

$\operatorname{max_{V}}$:

Let
$z\equiv[{\underline{z}},{\overline{z}}]\coloneqq\operatorname{max_{V}}(\hat{f},X)$.
As ${\overline{f}}$ is upper semicontinuous, then it must attain its maximum
over the compact set $X$ at some point $x_{0}\in X$. Thus, $\forall x\in
X:{\overline{f}}(x)\leq{\overline{f}}(x_{0})$. Suppose that $\epsilon>0$ is
given. By upper semicontinuity of ${\overline{f}}$, for each $x\in X$, there
exists a neighborhood $O_{x}$ such that $\forall y\in
O_{x}:{\overline{f}}(y)\leq{\overline{f}}(x)+\epsilon\leq{\overline{f}}(x_{0})+\epsilon$.
By referring to the sequence $(X_{i})_{i\in\mathbb{N}}$ defined in (23), and
using Scott continuity of $\hat{f}$, we deduce that, for each $x\in X$, there
exists a sufficiently large index $N_{x}\in\mathbb{N}$ that satisfies the
following property:

  * •

There exist $k_{x}\in\mathbb{N}$ and finitely many boxes
$\\{{B_{N_{x},j}}\mid{1\leq j\leq k_{x}}\\}$ such that:

$x\in{\left(\bigcup_{1\leq j\leq k_{x}}B_{N_{x},j}\right)}^{\circ}\text{ and
}\left(\bigcup_{1\leq j\leq k_{x}}B_{N_{x},j}\right)\subseteq O_{x},$ (24)

in which ${(\cdot)}^{\circ}$ denotes the interior operator, and:

$\forall
j\in\left\\{{1,\ldots,k_{x}}\right\\}:\quad{\overline{\hat{f}(B_{N_{x},j})}}\leq{\overline{f}}(x_{0})+\epsilon.$
(25)

The reason for considering $k_{x}$ boxes $\\{{B_{N_{x},j}}\mid{1\leq j\leq
k_{x}}\\}$—rather than just one box—is that, if at some iteration
$i_{0}\in\mathbb{N}$, the point $x$ lies on the boundary of a box in
$X_{i_{0}}$, then it will remain on the boundary of boxes in $X_{i}$ for all
$i\geq i_{0}$. In other words, it will not be in the interior of any boxes in
subsequent bisections. Of course, as each point $x\in X$ may be on the
boundaries of at most $2^{n}$ adjacent boxes in $X_{N_{x}}$, then we can
assume that $k_{x}\leq 2^{n}$.

For each $x\in X$, we define the neighborhood $U_{x}\coloneqq{(\bigcup_{1\leq
j\leq k_{x}}B_{N_{x},j})}^{\circ}$. As $X$ is compact, then it may be covered
by finitely many of such neighborhoods (say) $\\{{U_{x_{i}}}\mid{1\leq i\leq
m}\\}$. By taking $N=\max_{1\leq i\leq m}N_{x_{i}}$, and considering (25), we
deduce that, at the latest, at iteration $N$, we must have
${\overline{z_{i}}}\leq{\overline{z}}+\epsilon$.

One the other hand, as ${\underline{z}}\coloneqq\sup_{x\in
X}{\underline{f(x)}}$, for some $x^{\prime}_{0}\in X$ we must have
${\underline{f}}(x^{\prime}_{0})>{\underline{z}}-\epsilon/2$. Furthermore, as
${\underline{f}}$ is lower semicontinuous, then for some neighborhood
$X^{\prime}_{0}$ of $x^{\prime}_{0}$, we must have $\forall y\in
X^{\prime}_{0}:{\underline{f}}(y)>{\underline{f}}(x^{\prime}_{0})-\epsilon/2>{\underline{z}}-\epsilon$.
Similar to the previous argument, we can show that, after a sufficient number
of bisections, we must have ${\underline{z_{i}}}\geq{\underline{z}}-\epsilon$.

$\operatorname{max_{S}}$:

Assume that $x\in X\setminus\operatorname{max_{S}}(\hat{f},X)$. By Proposition
4.10, the set $\operatorname{max_{S}}(\hat{f},X)$ is compact. Hence, the set
$X\setminus\operatorname{max_{S}}(\hat{f},X)$ is relatively open in $X$ and
there exists a neighborhood $O_{x}$ of $x$ which satisfies $\forall y\in
O_{x}:{\overline{f}}(y)<{\underline{z}}$. Once again, we can show that, after
a sufficiently large number of bisections (say) $N_{x}$, for a finite number
of boxes, the relation (24) holds. This entails that all these boxes must be
discarded in line 9 of the algorithm after at most $N_{x}$ iterations.

∎

###  4.5 Computability

For computable analysis of the operators discussed in this work, we need an
effective structure on the underlying domains. To that end, we use the concept
of effectively given domains [75]. Specifically, we follow the approach taken
in [24, Section 3].

Assume that $(D,\sqsubseteq)$ is an $\omega$-continuous domain, with a
countable basis $B$ that is enumerated as follows:

$B=\left\\{{b_{0}=\bot,b_{1},\ldots,b_{n},\ldots}\right\\}.$ (26)

We say that the domain $D$ is _effectively given_ with respect to the
enumeration (26), if the set
$\\{{(i,j)\in\mathbb{N}\times\mathbb{N}}\mid{b_{i}\ll b_{j}}\\}$ is
recursively enumerable, in which $\ll$ is the way-below relation on $D$.

In the following proposition, computability is to be understood according to
TTE [84].

######  Proposition 4.14 (Computable elements and functions).

Let $D$ and $E$ be two effectively given domains, with enumerated bases
$B_{1}=\left\\{{d_{0},d_{1},\ldots,d_{n},\ldots}\right\\}$ and
$B_{2}=\left\\{{e_{0},e_{1},\ldots,e_{n},\ldots}\right\\}$, respectively. An
element $x\in D$ is computable iff the set $\\{{i\in\mathbb{N}}\mid{d_{i}\ll
x}\\}$ is recursively enumerable. A map $f:D\to E$ is computable iff
$\\{{(i,j)\in\mathbb{N}\times\mathbb{N}}\mid{e_{i}\ll f(d_{j})}\\}$ is
recursively enumerable.

###### Proof.

See [24, Proposition 3 and Theorem 9]. ∎

Arithmetic operators have computable domain-theoretic extensions, so do the
common activation functions (e. g., $\operatorname{\mathrm{ReLU}}$, $\tanh$,
or sigmoidal) and their Clarke-gradients [20, 33, 34, 21]. Our aim here is to
study the computability of $\operatorname{max_{V}}$ and
$\operatorname{max_{S}}$. From (17), we know that $\operatorname{max_{V}}$ and
$\operatorname{max_{S}}$ are not even monotone in their second argument,
while, by Theorem 4.11, if we fix $X\in\mathbb{K}\mathbb{R}^{n}_{\bot}$, then
the functionals
$\operatorname{max_{V}}(\cdot,X):({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
and
$\operatorname{max_{S}}(\cdot,X):({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to\mathbb{K}\mathbb{R}^{n}_{\bot}$
are Scott-continuous. By referring to Proposition 4.14, to study computability
of these two operators, we need to enumerate specific countable bases for the
domains
${\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$,
${\mathbb{I}}\mathbb{R}^{1}_{\bot}$, and $\mathbb{K}\mathbb{R}^{n}_{\bot}$. We
already know that $B_{{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$ and
$B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$ form countable bases for
${\mathbb{I}}\mathbb{R}^{1}_{\bot}$ and $\mathbb{K}\mathbb{R}^{n}_{\bot}$,
respectively, and it is straightforward to effectively enumerate them.

For the function space
${\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$,
we use the common step function construction. Recall the single-step functions
as defined in (15). Assume that $D$ and $E$ are two domains. For any $d\in D$
and $e\in E$, we must have:

$e\chi_{\twoheaduparrow{d}}(x)\coloneqq\left\\{\begin{array}[]{ll}e,&\text{if
}d\ll x,\\\ \bot,&\text{if }d\not\ll x.\end{array}\right.$

Note that, by Proposition 2.10, the set $\twoheaduparrow{d}$ is Scott-open. By
a step function we mean the supremum of a finite set of single-step functions.
If the domains $D$ and $E$ are bounded-complete, with countable bases
$B_{D}=\left\\{{d_{0},\ldots,d_{n},\ldots}\right\\}$ and
$B_{E}=\left\\{{e_{0},\ldots,e_{n},\ldots}\right\\}$, then the countable
collection:

$B_{D\Rightarrow E}\coloneqq\left\\{{\bigsqcup_{i\in
I}e_{i}\chi_{\twoheaduparrow{d_{i}}}}\;\vrule\;{I\text{ is
finite},\\{{e_{i}\chi_{\twoheaduparrow{d_{i}}}}\mid{i\in I}\\}\text{ is
bounded}}\right\\}$ (27)

forms a basis for the function space $D\Rightarrow E$ [1, Section 4].

Thus, we consider the bases $B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}$ and
$B_{{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$, and construct the countable basis
$B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$
for the function space
${\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}$,
according to (27). The key point is that, each step function in
$B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$
maps any given hyperbox with rational coordinates to a rational interval. We
formulate this as:

######  Proposition 4.15.

Each step function in
$B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$
maps any given element of $B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}}$ to an
element of $B_{{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$.

As such, we call the elements of
$B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$
rational step functions.

######  Theorem 4.16 (Computability).

Assume that $X\in B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$ is fixed, that is, $X$
is a finite union of hyperboxes with rational coordinates. Then, the
functionals
$\operatorname{max_{V}}(\cdot,X):({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to{\mathbb{I}}\mathbb{R}^{1}_{\bot}$
and
$\operatorname{max_{S}}(\cdot,X):({\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot})\to\mathbb{K}\mathbb{R}^{n}_{\bot}$
are computable.

###### Proof.

We prove computability according to Proposition 4.14. Let $f\in
B_{{\mathbb{I}}\mathbb{R}^{n}_{\bot}\Rightarrow{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$
be a rational step function.

$\operatorname{max_{V}}$:

Assume that $e\equiv[{\underline{e}},{\overline{e}}]\in
B_{{\mathbb{I}}\mathbb{R}^{1}_{\bot}}$, that is, $e$ is a rational interval.
To decide $e\ll\operatorname{max_{V}}(f,X)=[{\underline{z}},{\overline{z}}]$,
we run Algorithm 1. As $\twoheaduparrow{e}$ is Scott-open, then by
completeness (Theorem 4.13), after sufficiently large number of iterations
(say) $N_{f}$, we must have $\forall i\geq
N_{f}:e\ll[{{\underline{z_{i}}},{\overline{z_{i}}}}]$. By (3), the relation
$e\ll[{{\underline{z_{i}}},{\overline{z_{i}}}}]$ is equivalent to
${\underline{e}}<{\underline{z_{i}}}\leq{\overline{z_{i}}}<{\overline{e}}$. By
Proposition 4.15, both ${\underline{z_{i}}}$ and ${\overline{z_{i}}}$ are
rational numbers. Hence, deciding
$e\ll[{{\underline{z_{i}}},{\overline{z_{i}}}}]$ reduces to comparing rational
numbers, which is decidable. Thus, overall, the problem
$e\ll\operatorname{max_{V}}(f,X)$ is semi-decidable.

$\operatorname{max_{S}}$:

The proof for $\operatorname{max_{S}}$ is similar to that for
$\operatorname{max_{V}}$. Assume that $C\in
B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$. To decide
$C\ll\operatorname{max_{S}}(f,X)=K$, we run Algorithm 1. As
$\twoheaduparrow{C}$ is Scott-open, then by completeness (Theorem 4.13), after
sufficiently large number of iterations (say) $N_{f}$, we must have $\forall
i\geq N_{f}:C\ll\hat{B}_{i}$. By (3), the relation $C\ll\hat{B}_{i}$ is
equivalent to $\hat{B}_{i}\subseteq{C}^{\circ}$. As we have assumed that $C\in
B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$ and $X\in
B_{\mathbb{K}\mathbb{R}^{n}_{\bot}}$, then all the hyperboxes in $C$ and
$\hat{B}_{i}$ have rational coordinates. Therefore, deciding $C\ll\hat{B}_{i}$
reduces to comparing rational numbers, which is decidable. Thus, overall, the
problem $C\ll\operatorname{max_{S}}(f,X)$ is semi-decidable.

∎

###  4.6 Interval Arithmetic

We will use arbitrary-precision interval arithmetic for approximation of the
$\hat{L}$-derivative. Interval arithmetic [63] is performed over sets of
numbers, in contrast with the classical arithmetic which is performed over
individual numbers. This allows incorporating all the possible errors (e. g.,
round-off and truncations errors) into the result. As an example, note that:

$(x\in[{\underline{I}},{\overline{I}}])\wedge(y\in[{\underline{J}},{\overline{J}}])\implies
x+y\in[{\underline{I}}+{\underline{J}},{\overline{I}}+{\overline{J}}].$ (28)

In practice, inaccuracies enter into interval computations for a variety of
reasons. The following are the most relevant to the current article:

####  4.6.1 Finite Representations

In arbitrary-precision interval arithmetic, every real number $x$ is
represented as the limit of a sequence $(I_{i})_{i\in\mathbb{N}}$ of intervals
with dyadic end-points, i. e., $[x,x]=\bigcap_{i\in\mathbb{N}}I_{i}$. Assume
that we want to compute $\exp([1,1])=[\mathrm{e},\mathrm{e}]$, in which
$\mathrm{e}$ is the Euler number. As $\mathrm{e}$ is an irrational number, it
does not have a finite binary expansion. For instance, in MPFI [68]—the
arbitrary-precision interval library that we use—when the precision is set at
$10$, the result of $\exp([1,1])$ is given as $[2.718281828,2.718281829]$.

Furthermore, as the precision of the end-points of each interval must be set
first, inaccuracy is incurred even on simpler computations, such as (28). For
instance, if the bit size is set at $n$, then the result of
$[0,0]+[-2^{-2n},2^{-2n}]$ must be rounded outwards to $[-2^{-n},2^{-n}]$.

This source of inaccuracy is inevitable and, to a large extent, harmless. The
reason is that higher accuracy may be obtained by simply increasing the bit
size of the representation of the end-points.

####  4.6.2 Dependency

By default, interval computations do not take into account any dependencies
between parameters. For instance, assume that $x\in[0,1]$. By applying
interval subtraction we obtain $(x-x)\in[-1,1]$, which is a significant
overestimation as $x-x=0$, regardless of which interval $x$ belongs to.

Figure 2: Dependency problem: The network computes the constant function zero.
Hence, we must have $N^{\prime}(x)=0$ for all $x$. But the Clarke-gradient at
$x=0$ is overestimated as $[-1,1]$.

The dependency problem causes inaccuracies in estimations of the gradient as
well. As an example, consider the neural network of Figure 2. This network
computes the constant function $N(x)=0$, for which we must have $\forall
x\in\mathbb{R}:N^{\prime}(x)=0$. We estimate the Clarke-gradient of $N$ using
the chain rule, which, at $x=0$, leads to the following overestimation:

$\partial
N(0)\subseteq\partial\operatorname{\mathrm{ReLU}}(0)-\partial\operatorname{\mathrm{ReLU}}(0)=[0,1]-[0,1]=[-1,1].$

Our framework, however, is not affected by the dependency problem for general
position $\operatorname{\mathrm{ReLU}}$ networks (Theorem 4.19) or
differentiable networks (Theorem 4.21).

####  4.6.3 Wrapping Effect

Whenever a set which is not a hyperbox is overapproximated by a bounding
hyperbox, some accuracy is lost. This source of inaccuracy in interval
arithmetic is referred to as the wrapping effect [65]. As stated in
Proposition 4.2, the Clarke-gradient is, in general, a convex set. Hence, by
overapproximating the Clarke-gradient with a hyperbox, some accuracy will be
lost, as exemplified in Example 4.8.

Similar to the dependency problem, our framework is not affected by the
wrapping effect for general position $\operatorname{\mathrm{ReLU}}$ networks
(Theorem 4.19) and differentiable networks (Theorem 4.21).

###  4.7 General Position $\boldsymbol{\operatorname{\mathrm{ReLU}}}$ Networks

We have pointed out the trade-off between accuracy and efficiency in using the
hyperbox domain ${\mathbb{I}}\mathbb{R}^{n}_{\bot}$, and further sources of
inaccuracy in interval arithmetic. In this section, we show that, for a broad
class of $\operatorname{\mathrm{ReLU}}$ networks, no accuracy is lost by using
hyperboxes. Assume that $N:\mathbb{R}^{n}\to\mathbb{R}$ is a feedforward
$\operatorname{\mathrm{ReLU}}$ network with $k\geq 0$ hidden neurons (and an
arbitrary number of hidden layers). For convenience, we adopt the set-
theoretic definition of natural numbers [47] for which we have
$k=\left\\{{0,\ldots,k-1}\right\\}$. For each $i\in k$, let
$\hat{z}_{i}:\mathbb{R}^{n}\to\mathbb{R}$ be the pre-activation function of
the $i$-th hidden neuron and define:

$H_{i}\coloneqq\\{{x\in\mathbb{R}^{n}}\mid{\hat{z}_{i}(x)=0}\\}.$ (29)

The following definition, which we take from [49], is an adaptation of the
concept of hyperplane arrangements in general position [77]:

######  Definition 4.17.

A $\operatorname{\mathrm{ReLU}}$ network $N:\mathbb{R}^{n}\to\mathbb{R}$ with
$k\geq 0$ hidden neurons is said to be in _general position_ if, for every
subset $S\subseteq k$, the intersection $\bigcap_{i\in S}H_{i}$ is a finite
union of $(n-\lvert\thinspace{S}\thinspace\rvert)$-dimensional polytopes.

A $\operatorname{\mathrm{ReLU}}$ network $N:\mathbb{R}^{n}\to\mathbb{R}$ with
$k\geq 0$ hidden neurons and only one hidden layer is in general position if
for every $S\subseteq k$, the intersection $\bigcap_{i\in S}H_{i}$ is an
$(n-\lvert\thinspace{S}\thinspace\rvert)$-dimensional subspace of
$\mathbb{R}^{n}$. For instance, when $n=2$, the network is in general position
if no two hyperplanes $H_{i}$ and $H_{j}$ (with $i\neq j$) are parallel and no
three such hyperplanes meet at a point.

Recall that
$\tilde{L}(N):{\mathbb{I}}\mathbb{R}^{n}_{\bot}\to{\mathbb{I}}\mathbb{R}^{n}_{\bot}$
denotes the interval approximation of the $\hat{L}$-derivative of $N$ obtained
by applying the chain rule of Lemma 4.9. In general, we have:

$\forall x\in\mathbb{R}^{n}:\quad\partial
N(x)\subseteq\hat{L}(N)(x)\subseteq\tilde{L}(N)(\left\\{{x}\right\\}).$ (30)

Over general position networks, however, we obtain equality:

######  Lemma 4.18.

Let $N:\mathbb{R}^{n}\to\mathbb{R}$ be a general position
$\operatorname{\mathrm{ReLU}}$ network. Then, $\forall
x\in\mathbb{R}^{n}:\partial
N(x)=\hat{L}(N)(x)=\tilde{L}(N)(\left\\{{x}\right\\})$.

###### Proof.

The fact that $\forall x\in\mathbb{R}^{n}:\partial
N(x)=\tilde{L}(N)(\left\\{{x}\right\\})$ follows from [49, Theorem 2]. The
full claim now follows from (30). ∎

As a result, we have an effective procedure which, for any given
$\operatorname{\mathrm{ReLU}}$ network $N:[-M,M]^{n}\to\mathbb{R}$ in general
position, and any $X\in B_{\mathbb{K}[-M,M]^{n}}$, returns the Lipschitz
constant of $N$ over $X$, and the set of points in $X$ where the maximum
Lipschitz values are attained, to within any degree of accuracy:

######  Theorem 4.19.

Let $N:[-M,M]^{n}\to\mathbb{R}$ be a general position
$\operatorname{\mathrm{ReLU}}$ network and $X\in B_{\mathbb{K}[-M,M]^{n}}$.
Then:

$\left\\{\begin{array}[]{lcl}\operatorname{max_{V}}(\partial
N,X)&=&\operatorname{max_{V}}(\tilde{L}(N),X),\\\
\operatorname{max_{S}}(\partial
N,X)&=&\operatorname{max_{S}}(\tilde{L}(N),X),\end{array}\right.$

and both the maximum value and the maximum set can be obtained effectively to
any degree of accuracy.

###### Proof.

The result follows from Lemma 4.18, Theorem 4.13, and Theorem 4.16. ∎

The network of Figure 2 is not in general position, and we have shown that
Lemma 4.18 fails on that network. With respect to the Lebesgue measure over
the parameter space, however, almost every $\operatorname{\mathrm{ReLU}}$
network is in general position [49, Theorem 4]. In simple terms, _we have a
sound and effective procedure for Lipschitz analysis of networks, which
provides the results to within any given degree of accuracy on almost every
$\operatorname{\mathrm{ReLU}}$ network_.

###  4.8 Differentiable Networks

One of the main strengths of our domain-theoretic framework is that it can
handle differentiable and non-differentiable networks alike seamlessly. For
differentiable networks as well, we obtain results similar to those for
$\operatorname{\mathrm{ReLU}}$ networks in general position.

We know that when $N$ is classically differentiable at $x\in[-M,M]^{n}$, the
$L$-derivative and the classical derivative coincide [25]. Furthermore, for
classically differentiable functions, the chain rules of Proposition 4.3 and
Lemma 4.9 reduce to equalities which incur no overapproximation, because both
sides are single-valued. As a result, we obtain:

######  Lemma 4.20.

Let $N:[-M,M]^{n}\to\mathbb{R}$ be a differentiable feedforward network in
${\cal F}$. Then:

$\forall x\in[-M,M]^{n}:\quad\partial
N(x)=\hat{L}(N)(x)=\tilde{L}(N)(\left\\{{x}\right\\}).$

In other words, $\hat{L}(N)$ is an interval extension of $\partial N$.

######  Theorem 4.21.

Let $N:[-M,M]^{n}\to\mathbb{R}$ be a differentiable feedforward network in
${\cal F}$. Then:

$\left\\{\begin{array}[]{lcl}\operatorname{max_{V}}(\partial
N,X)&=&\operatorname{max_{V}}(\tilde{L}(N),X),\\\
\operatorname{max_{S}}(\partial
N,X)&=&\operatorname{max_{S}}(\tilde{L}(N),X),\end{array}\right.$

and both the maximum value and the maximum set can be obtained effectively to
any degree of accuracy.

##  5 Experiments

We have implemented our method for estimation of Lipschitz constant of neural
networks—which we call _lipDT_ —using two interval libraries: the arbitrary-
precision C++ library MPFI [68], and the fixed precision Python library
PyInterval.111https://github.com/taschini/pyinterval The latter is used for
fair time comparison with the other methods in Tables 1 and 3, which are all
in Python, i. e., ZLip [50], lipMIP [49], FastLip [86], SeqLip [82], and
CLEVER [85]. All the experiments have been conducted on a machine with an AMD
Ryzen 75800H 3.2GHz 8-Core processor and 16GB RAM under Ubuntu 20.04 LTS
environment. The source code is available on GitHub.222https://github.com/Can-
ZHOU/Theoretic-Robustness-NN

The main characteristic of our method is that it is fully validated. We begin
by demonstrating how our method can be used for verification of state of the
art methods from the literature through some experiments on very small
networks. If the state of the art methods turn out to be prone to inaccuracies
over small networks, then they most likely are prone to errors over large
networks as well.

Whenever a network has $m$ output neurons, we only consider the first output
neuron for analyses. We need to consider networks with more than one output to
avoid altering the methods from the literature which require more than one
output. For instance, in Experiment 1, we consider a randomly generated
$\operatorname{\mathrm{ReLU}}$ network of layer sizes $[2,2,2]$, with an input
domain straddling the intersection of the hyperplanes $H_{1}$ and $H_{2}$, as
defined in (29).

Table 1: Experiment 1: Random, $\operatorname{\mathrm{ReLU}}$, Layer sizes = [2, 2, 2]. Our method is called lipDT, which is the only case where both the method and the implementations are validated. Guarantee | Method | Estimate of Lip. Constant | Time (ms) | Rel. Err.
---|---|---|---|---
Fully validated | lipDT (C++) |  [0.32270569085905976, 0.32270569085905976] | 8.2890 | 0.00%
lipDT (Python) |  [0.32270569085905976, 0.32270569085905976] | 338.3580 | 0.00%
Method validated, implementation non-validated | ZLip | 0.041121020913124084 | 2.5075 | –87.26%
lipMIP | 0.041121020913124084 | 10.0348 | –87.26%
FastLip | 0.041121020913124084 | 0.8451 | –87.26%
Non-validated | SeqLip | 0.2129872590303421 | 1.9632 | –34.00%
CLEVER | 0.2815846800804138 | 1054.1022 | –12.74%

Assume that $z=[{\underline{z}},{\overline{z}}]$ is the interval returned by
lipDT. For each of the other methods, if $r$ is the returned estimate of the
Lipschitz constant, we calculate the relative error as follows:

$\text{relative
error}\coloneqq\left\\{\begin{array}[]{ll}\frac{r-{\underline{z}}}{{\underline{z}}},&\text{if
}r<{\underline{z}},\\\ 0,&\text{if }{\underline{z}}\leq
r\leq{\overline{z}},\\\ \frac{r-{\overline{z}}}{{\overline{z}}},&\text{if
}{\overline{z}}<r.\end{array}\right.$ (31)

As such, if the relative error is negative, then the value $r$ is definitely
an underapproximation of the true Lipschitz constant, and when the relative
error is positive, then $r$ is definitely an overapproximation of the true
Lipschitz constant. As reported in Tables 1 and 3, although lipMIP, ZLip, and
FastLip are based on validated methods, the implementations are not validated
due to floating-point errors. In Table 1, the estimates of the Lipschitz
constant returned by all the three methods are exactly the same, which
demonstrates that, for the input domain of Experiment 1, these methods all
follow the same computation path, which is subject to the same floating-point
errors.

For the randomly generated network of Experiment 1, the weights and biases are
shown in Figure 3. Based on these values, we calculate an input domain which
straddles the intersection of the hyperplanes $H_{1}$ and $H_{2}$. The center
and radius of the input domain are also given in Figure 3.

Weights $W_{1}$ from input neurons to the hidden layer:

$W_{1}=\begin{bmatrix}-0.5323450565338134766&0.2982044517993927002\\\
-0.4367777407169342041&0.1963265985250473022\end{bmatrix}$

Biases $b_{1}$ of the hidden layer:

$b_{1}=\begin{bmatrix}-0.3763356208801269531&-0.6647928357124328613\end{bmatrix}$

Weights $W_{2}$ from the hidden layer to the output layer:

$W_{2}=\begin{bmatrix}-0.3885447978973388672&0.4447681903839111328\\\
-0.3205857872962951660&0.2034754604101181030\end{bmatrix}$

Input domain:

$\text{center}=[x_{0},x_{1}]=[-4.832202221268014242,-7.364287590384273940],\quad\text{radius}=10^{-7}$

Figure 3: Weights, biases, and the input domain of the network of Experiment
1.

For lipDT to produce the results of Experiment 1, it had to perform five
iterations of the main loop of Algorithm 1. In each iteration, the most
important operation is that of the bisection, carried out in line 11 of
Algorithm 1. Table 2 shows the dynamics of the computation through successive
bisections.

Number of Bisections | Interval enclosure of the Lipschitz constant
---|---
0 | [0.0, 0.32270569085905976]
1 | [0.04112101957020187, 0.32270569085905976]
2 | [0.04112101957020187, 0.32270569085905976]
3 | [0.04112101957020187, 0.32270569085905976]
4 | [0.04112101957020187, 0.32270569085905976]
5 | [0.32270569085905976, 0.32270569085905976]
Table 2: The results of lipDT on Experiment 1 for the first 5 bisections.

The effects of the floating-point errors become less significant as the radius
of the input domain is increased, as demonstrated in Figure 4. We took the
same center as given in Figure 3, but changed the radius of the input domain
from $10^{-7}$ to $1$, and compared the results returned by all the methods.
ZLip, lipMIP, and FastLip give exactly the same results. Hence, we only draw
the results for ZLip, SeqLip, CLEVER, and lipDT. As can be seen, CLEVER and
ZLip—and as a consequence, the other methods except for SeqLip—all catch up
with lipDT eventually when the radius of the domain is large enough.

Figure 4: For very small input domains, lipDT is the only method which returns
the correct result.

###  5.1 Scalability: input dimension

When no bisections are required, lipDT can be useful for analysis of larger
networks as well, as shown in Experiment 2 on the MNIST dataset (Table 3). We
consider a $\operatorname{\mathrm{ReLU}}$ network of size $[784,10,10,2]$,
only focusing on the digits $1$ and $7$ in the output. The radius of the input
domain is $0.001$. For this experiment, the input domain is not situated on
the intersection of the hyperplanes, because lipDT would require bisecting
along 784 dimensions, resulting in $2^{784}$ new hyperboxes in each bisection,
which is prohibitively costly. Nonetheless, we observe that, even without any
bisection, lipDT maintains correctness, while ZLip, lipMIP, and FastLip still
underestimate the Lipschitz constant, albeit with small relative errors.

Table 3: Experiment 2: MNIST, $\operatorname{\mathrm{ReLU}}$, Layer sizes = [784, 10, 10, 2]. Guarantee | Method | Estimate of Lip. Constant | Time (ms) | Rel. Err.
---|---|---|---|---
Fully validated | lipDT C++ |  [0.34815343269413546, 0.34815343269413546] | 164.0400 | 0.00%
lipDT Python |  [0.34815343269412385, 0.3481534326941471] | 2388.6251 | 0.00%
Method validated, implementation non-validated | ZLip | 0.34815341234207153 | 7.1895 | –5.85e-06%
lipMIP | 0.3481534055299562 | 33.1788 | –7.80e-06%
FastLip | 0.34815341234207153 | 0.8844 | –5.85e-06%
Non-validated | SeqLip | 0.670467323694055 | 4.9305 | +92.58%
CLEVER | 0.3481535315513611 | 4308.8220 | +2.84e-05%

For a network with input dimension $n$, each bisection generates $2^{n}$ new
hyperboxes. As such, under favorable conditions, our method is applicable to
networks with input dimensions up to 10 or so. It is also possible to obtain
more accurate results on networks with higher input dimensions using
_influence analysis_ [83]. This can be achieved by a slight modification of
Algorithm 1. The influence analysis technique, however, is only applicable on
differentiable networks. Furthermore, this technique provides tighter bounds
on the Lipschitz constant in a reasonable time, but may still provide a loose
bound around the maximum set, i. e., the set where maximum Lipschitz constant
is attained. Another simple way of increasing efficiency is through
parallelization. The reason is that the algorithm is essentially performing a
search. Once several hyperboxes are generated, they can be searched in a
perfectly parallel way.

###  5.2 Scalability: depth and width

As we have mentioned previously, validated methods of Lipschitz estimation
(such as lipDT) cannot be scalable. Nonetheless, such methods are valuable for
purposes such as:

  1. (i)

verficiation of other Lipschitz estimation methods, as demonstrated in Table
1;

  2. (ii)

robustness analysis of low-dimensional networks, especially those that appear
in safety-critical applications, e. g., in power electronics [19], analog-to-
digital converters [81], and solution of PDE [69];

  3. (iii)

investigation of how Lipschitz constant is affected by factors such as under-
parametrization, over-parametrization, and through successive epochs of a
training process.

Item (iii) is particularly relevant for networks with low input dimensions,
but with depths and widths that can be very large relative to their input
dimensions. The reasons is that, although lipDT does not scale well with input
dimension of the network, it has a much better scalability with respect to the
depth and width of the network.

To investigate how lipDT scales with the number of layers (i. e., depth) of
the network, we generated some multilayer networks randomly with 4 input
neurons, 2 output neurons, and various number of hidden layers (from 1 to 10),
each containing 10 neurons. We took the input domain with center $[1,1,1,1]$
and radius $0.001$, and compared the timing of lipDT with the other methods,
except for CLEVER, which takes significantly longer than all the other
methods. As can be seen in Figure 5, lipDT scales quite well—almost
linearly—with the number of layers of the network.

Figure 5: lipDT scales (almost) linearly with respect to the number of layers.

Similarly, lipDT exhibits an almost linear scaling with respect to the width
of the network. We ran an experiment in which we considered networks with two
input neurons and one hidden layer of width $2^{n}$, where $n$ ranges from $1$
to $10$. The weights and biases were assigned randomly, and to make the
comparison fair, we performed the same number of bisections for each
experiment, making sure that the only contributing factor is the width. We
recorded the time, and the result is shown in Figure 6.

Figure 6: lipDT scales almost linearly with respect to the width of the
network.

###  5.3 Case study: benign overfitting

The traditional consensus in machine learning postulated that over-
parametrized models would have poor generalization properties. For instance,
in a regression task over a dataset of size $m$, if the capacity $K$ of the
model were significantly larger than $m$, then the model would interpolate the
dataset, and would oscillate erratically over points outside of the dataset.
In a fairly recent development, it has been observed that over-parametrized
models can indeed exhibit superior generalization compared with under-
parametrized ones, even over noisy data [64]. An extensive account of this
phenomenon—termed _‘benign overfitting’_ [5]—can be found in [6].

In this section, we investigate benign overfitting from the angle of Lipschitz
analysis. To that end, we consider a univariate regression task similar to
that depicted in [6, Figure 3.6]. The dataset consists of eleven points over
the interval $[-5,5]$. We train a set of $\operatorname{\mathrm{ReLU}}$
networks with one input neuron, one output neuron, and a hidden layer with
width ranging from $3$ to $3000$ (in steps of $3$). Even though lipDT scales
well with respect to the depth of the network, for this experiment, we
consider only networks with one hidden layer, and take the width of the single
hidden layer as the capacity of the network.

Each network is trained for $100\,000$ epochs. Then, the Lipschitz constant is
calculated over the interval $[-5,5]$ using lipDT. The result is shown in
Figure 7 (Left). To show the initial dynamics more clearly, the right hand
figure shows capacities only up to $200$. The curve exhibits similar
qualitative dynamics as that of the “double descent generalization curve” of
[6, Figure 3.5]. The under-parametrized regime extends from capacity $1$ up to
capacity $42$, during which the Lipschitz constant has an increasing and then
decreasing trend. Maximum Lipschitz constant is attained at capacity $21$.
Interpolation starts at capacity $42$, after which the Lipschitz constant
hovers around $4$ until it flattens after capacity $125$, and remains flat for
the remaining capacities up to $3000$.

Figure 7: Lipschitz constant versus capacity. Left: All the capacities up to
$3000$. Right: Capacities up to $200$. The under-parametrized regime extends
to capacity $42$, where interpolation starts, and after which, the curve
remains almost flat. Maximum Lipschitz constant is attained at capacity $21$.

Figure 8 depicts the output of neural networks with capacities $3$ and $21$
(under-parametrized), $42$ (interpolation threshold) and $200$ (over-
parametrized). The slope of the linear segment connecting the third and fourth
data points (from left) is roughly $4$. By the mean value theorem, it is clear
that once interpolation starts, the Lipschitz constant may not go below this
slope. As can be seen in Figure 7, in the over-parametrized regime, the
Lipschitz constant stays very close to $4$.

We also investigate how the Lipschitz constant is affected through training.
The results are shown in Figure 9. In particular, when the capacity is farther
away from the interpolation threshold—e. g., for capacities $3$ and $200$—the
Lipschitz curve flattens at earlier epochs. In contrast, for capacities $21$
and $42$, optimization seems a slower and more demanding process, as the curve
begins to flatten at much later epochs. We point out, however, that the
optimization algorithm is randomized, and the network is initialized randomly.
Hence, the exact dynamics of the Lipschitz constant versus epoch number—in
particular, at what stage the curve begins to flatten—are subject to the
initial weights and biases and the optimization process.

Figure 8: The dashed black curve is a piecewise linear interpolation of the
data points. The blue curves show the output of the neural networks.
Interpolation starts at capacity $42$. The highest Lipschitz constant is
attained at capacity $21$. Also, see, [6, Figure 3.6].

Figure 9: Dynamics of the Lipschitz constant through successive epochs. The
network with three neurons reaches its best loss much earlier than $100\,000$
epochs, so we stop the training at $20000$ epochs.

####  5.3.1 Validated versus non-validated Lipschitz analysis

We ran our experiments using the validated method lipDT, and the non-validated
statistical method CLEVER. The Lipschitz estimates returned by CLEVER are
almost always a slight under-estimate of those returned by lipDT. To be more
precise, the average relative error of CLEVER is around $-1\%$, calculated
according to (31). In most cases, qualitatively, the results of CLEVER and
lipDT are similar. Some outliers, however, were observed. For instance, Figure
10 (Left) shows the graph of the result returned by CLEVER when measuring the
Lipschitz constant of the $\operatorname{\mathrm{ReLU}}$ network with $3000$
neurons, up to epoch $2500$. At epoch $251$, CLEVER estimated the Lipschitz
constant to be around $127$, while for all the other epochs, the Lipschitz
constant was estimated below $4$.

Figure 10: Left: CLEVER returns an outlier at epoch $251$. Right: Comparison
with lipDT (for the first $350$ epochs) reveals that the Lipschitz constant
does not have a jump.

Without recourse to a validated method, it is not possible to determine
whether the Lipschitz constant truly has a jump at epoch $251$ or not. The
reason is that the commonly used optimization methods are randomized, which
makes it possible to obtain genuine outliers at certain epochs. Nonetheless,
running the same experiment with lipDT reveals that, at epoch $251$, CLEVER
indeed provides the wrong result. The interval returned by lipDT at epoch
$251$ is:

$[2.445572934050687,2.445572934050933].$

Figure 10 (Right) depicts the results (for the first $350$ epochs) of the
comparison between lipDT and CLEVER.

While in experiments of Table 1 and Figure 4, the input domains were chosen
deliberately to be small and to straddle the intersection of hyperplanes, in
the experiment of Figure 10, the search was performed over the entire interval
of $[-5,5]$. This demonstrates that non-validated statistical methods may
return misleading results even when the samples are taken from a large and
seemingly non-critical part of a very low-dimensional input space.

##  6 Concluding remarks

Continuous domains are powerful models of computation and provide a semantic
model for the concept of approximation. We have demonstrated this claim by
presenting a framework for robustness analysis of neural networks. The main
focus in the current article has been laying the theoretical foundation and
obtaining fundamental results on soundness, completeness, and computability.
We also demonstrated how straightforward it is to translate our algorithms
into fully validated implementations using interval arithmetic. In particular,
we devised an algorithm for validated computation of Lipschitz constants of
feedforward networks, which is sound and complete. Based on the foundation
laid in this paper, our next step will be extending the framework to RNN.

The maximization algorithm of Algorithm 1 has an imperative style, which is in
line with the common practice in machine learning. Continuous domains,
however, are best known as models of lambda calculus and it is common to adopt
a functional style in language design for real number computation [33, 34,
17]. Thus, another direction for future work is the design of a functional
algorithm for maximization which suits our non-algebraic domain framework.

## References

  * [1] S. Abramsky and A. Jung  “Domain Theory”  In _Handbook of Logic in Computer Science_ 3 Clarendon Press, Oxford, 1994, pp. 1–168
  * [2] Samson Abramsky  “Abstract Interpretation, Logical Relations, and Kan Extensions”  In _Journal of Logic and Computation_ 1.1, 1990, pp. 5–40
  * [3] Fernando Albiac and Nigel J. Kalton  “Topics in Banach Space Theory”  Springer, 2006
  * [4] Alexandre Araujo, Benjamin Négrevergne, Yann Chevaleyre and Jamal Atif  “On Lipschitz Regularization of Convolutional Layers using Toeplitz Matrix Theory”  In _AAAI_ , 2021
  * [5] Peter L. Bartlett, Philip M. Long, Gábor Lugosi and Alexander Tsigler  “Benign overfitting in linear regression”  In _Proceedings of the National Academy of Sciences_ 117.48, 2020, pp. 30063–30070  DOI: 10.1073/pnas.1907378117
  * [6] Mikhail Belkin  “Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation”  In _Acta Numerica_ 30 Cambridge University Press, 2021, pp. 203–248  DOI: 10.1017/S0962492921000039
  * [7] Aritra Bhowmick, Meenakshi D’Souza and G. Srinivasa Raghavan  “LipBaB: Computing Exact Lipschitz Constant of ReLU Networks”  In _Artificial Neural Networks and Machine Learning – ICANN 2021_ Cham: Springer International Publishing, 2021, pp. 151–162  DOI: 10.1007/978-3-030-86380-7˙13
  * [8] Geoffrey L. Burn, Chris Hankin and Samson Abramsky  “Strictness analysis for higher-order functions”  In _Science of Computer Programming_ 7, 1986, pp. 249–278
  * [9] Nicholas Carlini and David A. Wagner  “Towards Evaluating the Robustness of Neural Networks”  In _2017 IEEE Symposium on Security and Privacy, SP 2017, San Jose, CA, USA, May 22-26, 2017_ IEEE Computer Society, 2017, pp. 39–57  DOI: 10.1109/SP.2017.49
  * [10] Swarat Chaudhuri, Sumit Gulwani and Roberto Lublinerman  “Continuity and Robustness of Programs”  In _Commun. ACM_ 55.8, 2012, pp. 107–115
  * [11] Swarat Chaudhuri, Sumit Gulwani, Roberto Lublinerman and Sara Navidpour  “Proving Programs Robust”  In _Proceedings of the 19th ACM SIGSOFT Symposium and the 13th European Conference on Foundations of Software Engineering_ , ESEC/FSE ’11  Association for Computing Machinery, 2011, pp. 102–112
  * [12] Tong Chen, Jean B. Lasserre, Victor Magron and Edouard Pauwels  “Semialgebraic Optimization for Lipschitz Constants of ReLU Networks”  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6–12, 2020, virtual_ , 2020
  * [13] F. H. Clarke, Yu .S. Ledyaev, R. J. Stern and R. R. Wolenski  “Nonsmooth Analysis and Control Theory”  Springer, 1998
  * [14] Frank H. Clarke  “Optimization and Nonsmooth Analysis”, Classics in Applied Mathematics  Society for IndustrialApplied Mathematics, 1990
  * [15] Patrick Cousot and Radhia Cousot  “Abstract Interpretation: A Unified Lattice Model for Static Analysis of Programs by Construction or Approximation of Fixpoints”  In _Proceedings of the 4th ACM SIGACT-SIGPLAN Symposium on Principles of Programming Languages_ , POPL ’77, 1977, pp. 238–252
  * [16] Pietro Di Gianantonio  “Real Number Computability and Domain Theory”  In _Information and Computation_ 127.1, 1996, pp. 11–25
  * [17] Pietro Di Gianantonio and Abbas Edalat  “A Language for Differentiable Functions”  In _Foundations of Software Science and Computation Structures_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2013, pp. 337–352  DOI: 10.1007/978-3-642-37075-5˙22
  * [18] Thomas G. Dietterich  “Steps Toward Robust Artificial Intelligence”  In _AI Magazine_ 38.3, 2017, pp. 3–24  DOI: 10.1609/aimag.v38i3.2756
  * [19] Tomislav Dragičević, Patrick Wheeler and Frede Blaabjerg  “Artificial Intelligence Aided Automated Design for Reliability of Power Electronic Systems”  In _IEEE Transactions on Power Electronics_ 34.8, 2019, pp. 7161–7171  DOI: 10.1109/TPEL.2018.2883947