11institutetext: Stony Brook University, Stony Brook, NY, USA, 11email:
<EMAIL_ADDRESS>22institutetext: Aalborg
University, Aalborg, Denmark, 22email<EMAIL_ADDRESS>33institutetext:
Technichal University of Munich, Garching, Germany, 33email<EMAIL_ADDRESS>
# Open- and Closed-Loop Neural Network Verification using Polynomial Zonotopes

Niklas Kochdumper 11 Christian Schilling 22 Matthias Althoff 33 and Stanley
Bak 11

###### Abstract

We present a novel approach to efficiently compute tight non-convex enclosures
of the image through neural networks with ReLU, sigmoid, or hyperbolic tangent
activation functions. In particular, we abstract the input-output relation of
each neuron by a polynomial approximation, which is evaluated in a set-based
manner using polynomial zonotopes. While our approach can also can be
beneficial for open-loop neural network verification, our main application is
reachability analysis of neural network controlled systems, where polynomial
zonotopes are able to capture the non-convexity caused by the neural network
as well as the system dynamics. This results in a superior performance
compared to other methods, as we demonstrate on various benchmarks.

###### Keywords:

Neural network verification neural network controlled systems reachability
analysis polynomial zonotopes formal verification.

##  1 Introduction

While previously artificial intelligence was mainly used for soft applications
such as movie recommendations [9], facial recognition [23], or chess computers
[11], it is now also increasingly applied in safety-critical applications,
such as autonomous driving [32], human-robot collaboration [27], or power
system control [5]. In contrast to soft applications, where failures usually
only have minor consequences, failures in safety-critical applications in the
worst case result in loss of human lives. Consequently, in order to prevent
those failures, there is an urgent need for efficient methods that can verify
that the neural networks used for artificial intelligence function correctly.
Verification problems involving neural networks can be grouped into two main
categories:

  * •

Open-loop verification: Here the task is to check if the output of the neural
network for a given input set satisfies certain properties. With this setup
one can for example prove that a neural network used for image classification
is robust against a certain amount of noise on the image.

  * •

Closed-loop verification: In this case the neural network is used as a
controller for a dynamical system, e.g., to steer the system to a given goal
set while avoiding unsafe regions. The safety of the controlled system can be
verified using reachability analysis.

For both of the above verification problems, the most challenging step is to
compute a tight enclosure of the image through the neural network for a given
input set. Due to the high expressiveness of neural networks, their images
usually have complex shapes, so that convex enclosures are often too
conservative for verification. In this work, we show how to overcome this
limitation with our novel approach for computing tight non-convex enclosures
of images through neural networks using polynomial zonotopes.

###  1.1 State of the Art

We first summarize the state of the art for open-loop neural network
verification followed by reachability analysis for neural network controlled
systems. Many different set representations have been proposed for computing
enclosures of the image through a neural network, including intervals [43],
polytopes [38], zonotopes [34], star sets [40], and Taylor models [21]. For
neural networks with ReLU activation functions, it is possible to compute the
exact image. This can be either achieved by recursively partitioning the input
set into piecewise affine regions [42], or by propagating the initial set
through the network using polytopes [38, 48] or star sets [40], where the set
is split at all neurons that are both active or inactive. In either case the
exact image is in the worst case given as a union of $2^{v}$ convex sets, with
$v$ being the number of neurons in the network. To avoid this high
computational complexity for exact image computation, most approaches compute
a tight enclosure instead using an abstraction of the neural network. For ReLU
activation functions one commonly used abstraction is the triangle relaxation
[15] (see Fig. 1), which can be conveniently integrated into set propagation
using star sets [40]. Another possibility is to abstract the input-output
relation by a zonotope (see Fig. 1), which is possible for ReLU, sigmoid, and
hyperbolic tangent activation functions [34]. One can also apply Taylor model
arithmetic [26] to compute the image through networks with sigmoid and
hyperbolic tangent activation [21], which corresponds to an abstraction of the
input-output relation by a Taylor series expansion. In order to better capture
dependencies between different neurons, some approaches also abstract the
input-output relation of multiple neurons at once [36, 28].

Figure 1: Triangle relaxation (left), zonotope abstraction (middle), and
polynomial zonotope abstraction (right) of the ReLU activation function.

While computation of the exact image is infeasible for large networks, the
enclosures obtained by abstractions are often too conservative for
verification. To obtain complete verifiers, many approaches therefore use
branch and bound strategies [7] that split the input set and/or single neurons
until the specification can either be proven or a counterexample is found. For
computational reasons branch and bound strategies are usually combined with
approaches that are able to compute rough interval bounds for the neural
network output very fast. Those bounds can for example be obtained using
symbolic intervals [43] that store linear constraints on the variables in
addition to the interval bounds to preserve dependencies. The DeepPoly
approach [35] uses a similar concept, but applies a back-substitution scheme
to obtain tighter bounds. With the FastLin method [45] linear bounds for the
overall network can be computed from linear bounds for the single neurons. The
CROWN approach [49] extends this concept to linear bounds with different
slopes as well as quadratic bounds. Several additional improvements for the
CROWN approach have been proposed, including slope optimization using gradient
descent [47] and efficient ReLU splitting [44]. Instead of explicitly
computing the image, many approaches also aim to verify the specification
directly using SMT solvers [30, 22], mixed-integer linear programming [8, 37],
semidefinite programming [31], and convex optimization [24].

For reachability analysis of neural network controlled systems one has to
compute the set of control inputs in each control cycle, which is the image of
the current reachable set through the neural network controller. Early
approaches compute the image for ReLU networks exactly using polytopes [46] or
star sets [39]. Since in this case the number of coexisting sets grows rapidly
over time, these approaches have to unite sets using convex hulls [46] or
interval enclosures [39], which often results in large over-approximations. If
template polyhedra are used as a set representation, reachability analysis for
neural network controlled systems with discrete-time plants reduces to the
task of computing the maximum output along the template directions [12], which
can be done efficiently. Neural network controllers with sigmoid and
hyperbolic tangent activation functions can be converted to an equivalent
hybrid automaton [20], which can be combined with the dynamics of the plant
using the automaton product. However, since each neuron is represented by an
additional state, the resulting hybrid automaton is very high-dimensional,
which makes reachability analysis challenging. Some approaches approximate the
overall network with a polynomial function [14, 18] using polynomial
regression based on samples [14] and Bernstein polynomials [18]. Yet another
class of methods [10, 21, 33, 41] employs abstractions of the input-output
relation for the neurons to compute the set of control inputs using intervals
[10], star sets [41], Taylor models[21], and a combination of zonotopes and
Taylor models[33]. Common tools for reachability analysis of neural network
controlled systems are JuliaReach [6], NNV [41], POLAR [19], ReachNN* [16],
RINO [17], Sherlock [13], Verisig [20], and Verisig 2.0 [21], where JuliaReach
uses zonotopes for neural network abstraction [33], NVV supports multiple set
representations, ReachNN* applies the Bernstein polynomial method [18], POLAR
approximates single neurons by Bernstein polynomials [19], RINO computes
interval inner- and outer-approximations [17], Sherlock uses the polynomial
regression approach [14], Verisig performs the conversion to a hybrid
automaton [20], and Verisig 2.0 uses the Taylor model based neural network
abstraction method [21].

###  1.2 Overview

In this work, we present a novel approach for computing tight non-convex
enclosures of images through neural networks with ReLU, sigmoid, or hyperbolic
tangent activation functions. The high-level idea is to approximate the input-
output relation of each neuron by a polynomial function, which results in the
abstraction visualized in Fig. 1. Since polynomial zonotopes are closed under
polynomial maps, the image through this function can be computed exactly,
yielding a tight enclosure of the image through the overall neural network.
The remainder of this paper is structured as follows: After introducing some
preliminaries in Sec. 2, we present our approach for computing tight
enclosures of images through neural networks in Sec. 3. Next, we show how to
utilize this result for reachability analysis of neural network controlled
systems in Sec. 4. Afterwards, in Sec. 5, we introduce some special operations
on polynomial zonotopes that we require for image and reachable set
computation, before we finally demonstrate the performance of our approach on
numerical examples in Sec. 6.

###  1.3 Notation

Sets are denoted by calligraphic letters, matrices by uppercase letters, and
vectors by lowercase letters. Given a vector $b\in\mathbb{R}^{n}$, $b_{(i)}$
refers to the $i$-th entry. Given a matrix $A\in\mathbb{R}^{o\times n}$,
$A_{(i,\cdot)}$ represents the $i$-th matrix row, $A_{(\cdot,j)}$ the $j$-th
column, and $A_{(i,j)}$ the $j$-th entry of matrix row $i$. The concatenation
of two matrices $C$ and $D$ is denoted by $[C~{}D]$, and
$I_{n}\in\mathbb{R}^{n\times n}$ is the identity matrix. The symbols
$\mathbf{0}$ and $\mathbf{1}$ represent matrices of zeros and ones of proper
dimension, the empty matrix is denoted by $[~{}]$, and $\mathrm{diag}(a)$
returns a diagonal matrix with $a\in\mathbb{R}^{n}$ on the diagonal. Given a
function $f(x)$ defined as $f:~{}\mathbb{R}\to\mathbb{R}$, $f^{\prime}(x)$ and
$f^{\prime\prime}(x)$ denote the first and second derivative with respect to
$x$. The left multiplication of a matrix $A\in\mathbb{R}^{o\times n}$ with a
set $\mathcal{S}\subset\mathbb{R}^{n}$ is defined as
$A\,\mathcal{S}:=\\{A\,s~{}|~{}s\in\mathcal{S}\\}$, the Minkowski addition of
two sets $\mathcal{S}_{1}\subset\mathbb{R}^{n}$ and
$\mathcal{S}_{2}\subset\mathbb{R}^{n}$ is defined as
$\mathcal{S}_{1}\oplus\mathcal{S}_{2}:=\\{s_{1}+s_{2}~{}|~{}s_{1}\in\mathcal{S}_{1},s_{2}\in\mathcal{S}_{2}\\}$,
and the Cartesian product of two sets $\mathcal{S}_{1}\subset\mathbb{R}^{n}$
and $\mathcal{S}_{2}\subset\mathbb{R}^{o}$ is defined as
$\mathcal{S}_{1}\times\mathcal{S}_{2}:=\big{\\{}[s_{1}^{T}~{}s_{2}^{T}]^{T}~{}|~{}s_{1}\in\mathcal{S}_{1},s_{2}\in\mathcal{S}_{2}\big{\\}}$.
We further introduce an $n$-dimensional interval as
$\mathcal{I}:=[l,u],~{}\forall i~{}l_{(i)}\leq
u_{(i)},~{}l,u\in\mathbb{R}^{n}$.

##  2 Preliminaries

Let us first introduce some preliminaries required throughout the paper. While
the concepts presented in this work can equally be applied to more complex
network architectures, we focus on feed-forward neural networks for
simplicity:

###### Definition 1

(Feed-forward neural network) A feed-forward neural network with $\kappa$
hidden layers consists of weight matrices $W_{i}\in\mathbb{R}^{v_{i}\times
v_{i-1}}$ and bias vectors $b_{i}\in\mathbb{R}^{v_{i}}$ with
$i\in\\{1,\dots,\kappa+1\\}$ and $v_{i}$ denoting the number of neurons in
layer $i$. The output $y\in\mathbb{R}^{v_{\kappa+1}}$ of the neural network
for the input $x\in\mathbb{R}^{v_{0}}$ is

$y:=y_{\kappa+1}~{}~{}\text{with}~{}~{}y_{0}=x,~{}~{}y_{i(j)}=\mu\bigg{(}\sum_{k=1}^{v_{i-1}}W_{i(j,k)}\,y_{i-1(k)}+b_{i(j)}\bigg{)},~{}i=1,\dots,\kappa+1,$

where $\mu:~{}\mathbb{R}\to\mathbb{R}$ is the activation function.

In this paper we consider ReLU activations $\mu(x)=\max(0,x)$, sigmoid
activations $\mu(x)=\sigma(x)=1/(1+e^{-x})$, and hyperbolic tangent
activations $\mu(x)=\tanh(x)=(e^{x}-e^{-x})/(e^{x}+e^{-x})$. Moreover, neural
networks often do not apply activation functions on the output neurons, which
corresponds to using the identity map $\mu(x)=x$ for the last layer. The image
$\mathcal{Y}$ through a neural network is defined as the set of outputs for a
given set of inputs $\mathcal{X}_{0}$, which is according to Def. 1 given as

$\mathcal{Y}=\bigg{\\{}y_{\kappa+1}~{}\bigg{|}~{}y_{0}\in\mathcal{X}_{0},~{}\forall
i\in\\{1,\dots,\kappa+1\\}:~{}y_{i(j)}=\mu\bigg{(}\sum_{k=1}^{v_{i-1}}W_{i(j,k)}\,y_{i-1(k)}+b_{i(j)}\bigg{)}\bigg{\\}}.$

We present a novel approach for tightly enclosing the image through a neural
network by a polynomial zonotope [2], where we use the sparse representation
of polynomial zonotopes [25]-2-2-2In contrast to [25, Def. 1], we explicitly
do not integrate the constant offset $c$ in $G$. Moreover, we omit the
identifier vector used in [25] for simplicity:

Figure 2: Step-by-step construction of the polynomial zonotope from Example 1.

###### Definition 2

(Polynomial zonotope) Given a constant offset $c\in\mathbb{R}^{n}$, a
generator matrix of dependent generators $G\in\mathbb{R}^{n\times h}$, a
generator matrix of independent generators $G_{I}\in\mathbb{R}^{n\times q}$,
and an exponent matrix $E\in\mathbb{N}_{0}^{p\times h}$, a polynomial zonotope
$\mathcal{PZ}\subset\mathbb{R}^{n}$ is defined as

$\mathcal{PZ}:=\bigg{\\{}c+\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\bigg{)}G_{(\cdot,i)}+\sum_{j=1}^{q}\beta_{j}G_{I(\cdot,j)}~{}\bigg{|}~{}\alpha_{k},\beta_{j}\in[-1,1]\bigg{\\}}.$

The scalars $\alpha_{k}$ are called  dependent factors since a change in their
value affects multiplication with multiple generators. Analogously, the
scalars $\beta_{j}$ are called independent factors because they only affect
the multiplication with one generator. For a concise notation we use the
shorthand $\mathcal{PZ}=\langle c,G,G_{I},E\rangle_{PZ}$.

Let us demonstrate polynomial zonotopes by an example:

###### Example 1

The polynomial zonotope

$\mathcal{PZ}=\left\langle\begin{bmatrix}4\\\
4\end{bmatrix},\begin{bmatrix}2&~{}1~{}&2\\\
0&~{}2~{}&2\end{bmatrix},\begin{bmatrix}1\\\
0\end{bmatrix},\begin{bmatrix}1&~{}0~{}&3\\\
0&1&1\end{bmatrix}\right\rangle_{PZ}$

defines the set

$\mathcal{PZ}=\bigg{\\{}\begin{bmatrix}4\\\ 4\end{bmatrix}+\begin{bmatrix}2\\\
0\end{bmatrix}\alpha_{1}+\begin{bmatrix}1\\\
2\end{bmatrix}\alpha_{2}+\begin{bmatrix}2\\\
2\end{bmatrix}\alpha_{1}^{3}\alpha_{2}+\begin{bmatrix}1\\\
0\end{bmatrix}\beta_{1}~{}\bigg{|}~{}\alpha_{1},\alpha_{2},\beta_{1}\in[-1,1]\bigg{\\}}.$

The construction of this polynomial zonotope is visualized in Fig. 2.

##  3 Image Enclosure

We now present our novel approach for computing tight non-convex enclosures of
images through neural networks. The general concept is to approximate the
input-output relation of each neuron by a polynomial function, the image
through which can be computed exactly since polynomial zonotopes are closed
under polynomial maps. For simplicity, we focus on quadratic approximations
here, but the extension to polynomials of higher degree is straightforward.

The overall procedure for computing the image is summarized in Alg. 1, where
the computation proceeds layer by layer. For each neuron in the current layer
$i$ we first calculate the corresponding input set in Line 5. Next, in Line 6,
we compute a lower and an upper bound for the input to the neuron. Using these
bounds we then calculate a quadratic approximation for the neuron’s input-
output relation in Line 7. This approximation is evaluated in a set-based
manner in Line 8. The resulting polynomial zonotope $\langle
c_{q},G_{q},G_{I,q},E_{q}\rangle_{PZ}$ forms the $j$-th dimension of the set
$\mathcal{PZ}$ representing the output of the whole layer (see Line 9 and Line
12). To obtain a formally correct enclosure, we have to account for the error
made by the approximation. We therefore compute the difference between the
activation function and the quadratic approximation in Line 10 and add the
result to the output set in Line 12. By repeating this procedure for all
layers, we finally obtain a tight enclosure of the image through the neural
network. A demonstrating example for Alg. 1 is shown in Fig. 3.

For ReLU activations the quadratic approximation only needs to be calculated
if $l<0\wedge u>0$ since we can use the exact input-output relations $g(x)=x$
and $g(x)=0$ if $l\geq 0$ or $u\leq 0$ holds. Due to the evaluation of the
quadratic map defined by $g(x)$, the representation size of the polynomial
zonotope $\mathcal{PZ}$ increases in each layer. For deep neural networks it
is therefore advisable to repeatedly reduce the representation size after each
layer using order reduction [25, Prop. 16]. Moreover, one can also apply the
compact operation described in [25, Prop. 2] after each layer to remove
potential redundancies from $\mathcal{PZ}$. Next, we explain the approximation
of the input-output relation as well as the computation of the approximation
error in detail.

Algorithm 1 Enclosure of the image through a neural network

Require: Neural network with weight matrices $W_{i}$ and bias vectors $b_{i}$,
initial set $\mathcal{X}_{0}$.

Ensure: Tight enclosure $\mathcal{PZ}\supseteq\mathcal{Y}$ of the image
$\mathcal{Y}$.

1:$\mathcal{PZ}\leftarrow\mathcal{X}_{0}$

2:for $i\leftarrow 1$ to $\kappa+1$ do (loop over all layers)

3:
$c\leftarrow\mathbf{0},~{}G\leftarrow\mathbf{0},~{}G_{I}\leftarrow\mathbf{0},~{}\underline{d}\leftarrow\mathbf{0},~{}\overline{d}\leftarrow\mathbf{0}$

4: for $j\leftarrow 1$ to $v_{i}$ do (loop over all neurons in the layer)

5: $\mathcal{PZ}_{j}\leftarrow W_{i(j,\cdot)}\mathcal{PZ}\oplus b_{i(j)}$ (map
with weight matrix and bias using (5))

6: $l,u\leftarrow$ lower and upper bound for $\mathcal{PZ}_{j}$ according to
Prop. 1

7: $g(x)=a_{1}\,x^{2}+a_{2}\,x+a_{3}\leftarrow$ quad. approx. on $[l,u]$
according to Sec. 3.1

8: $\langle c_{q},G_{q},G_{I,q},E_{q}\rangle_{PZ}\leftarrow$ image of
$\mathcal{PZ}_{j}$ through $g(x)$ according to Prop. 2

9: $c_{(j)}\leftarrow c_{q},~{}G_{(j,\cdot)}\leftarrow
G_{q},~{}G_{I(j,\cdot)}\leftarrow G_{I,q},~{}E\leftarrow E_{q}$ (add to output
set)

10: $\underline{d}_{(j)},\overline{d}_{(j)}\leftarrow$ difference between
$g(x)$ and activation function acc. to Sec. 3.2

11: end for

12: $\mathcal{PZ}\leftarrow\langle
c,G,G_{I},E\rangle_{PZ}\oplus[\underline{d},\overline{d}]$ (add approximation
error using (6))

13:end for

###  3.1 Activation Function Approximation

The centerpiece of our algorithm for computing the image of a neural network
is the approximation of the input-output relation defined by the activation
function $\mu(x)$ with a quadratic expression
$g(x)=a_{1}\,x^{2}+a_{2}\,x+a_{3}$ (see Line 7 of Alg. 1). In this section we
present multiple possibilities to obtain good approximations.

#### Polynomial Regression

For polynomial regression we uniformly select $N$ samples $x_{i}$ from the
interval $[l,u]$ and then determine the polynomial coefficients
$a_{1},a_{2},a_{3}$ by minimizing the average squared distance between the
activation function and the quadratic approximation:
$\min_{a_{1},a_{2},a_{3}}\frac{1}{N}\sum_{i=1}^{N}\big{(}\mu(x_{i})-a_{1}\,x_{i}^{2}-a_{2}\,x_{i}-a_{3}\big{)}^{2}.$
(1) It is well known that the optimal solution to (1) is
$\begin{bmatrix}a_{1}\\\ a_{2}\\\
a_{3}\end{bmatrix}=A^{\dagger}b~{}~{}\text{with}~{}~{}A=\begin{bmatrix}x_{1}^{2}&x_{1}&1\\\
\vdots&\vdots&\vdots\\\
x_{N}^{2}&x_{N}&1\end{bmatrix},~{}~{}b=\begin{bmatrix}\mu(x_{1})\\\ \vdots\\\
\mu(x_{N})\end{bmatrix},$ where $A^{\dagger}=(A^{T}A)^{-1}A^{T}$ is the Moore-
Penrose inverse of matrix A. For the numerical experiments in this paper we
use $N=10$ samples.

#### Closed-Form Expression

For ReLU activations a closed-form expression for a quadratic approximation
can be obtained by enforcing the conditions $g(l)=0$, $g^{\prime}(l)=0$, and
$g(u)=u$. The solution to the corresponding equation system
$a_{1}\,l^{2}+a_{2}\,l+a_{3}=0$, $2a_{1}l+a_{2}=0$,
$a_{1}\,u^{2}+a_{2}\,u+a_{3}=u$ is
$a_{1}=\frac{u}{(u-l)^{2}},~{}~{}a_{2}=\frac{-2lu}{(u-l)^{2}},~{}~{}a_{3}=\frac{u^{2}(2l-u)}{(u-l)^{2}}+u,$
which results in the enclosure visualized in Fig. 1. This closed-form
expression is very precise if the interval $[l,u]$ is close to being symmetric
with respect to the origin ($|l|\approx|u|$), but becomes less accurate if one
bound is significantly larger than the other ($|u|\gg|l|$ or $|l|\gg|u|$).

#### Taylor Series Expansion

For sigmoid and hyperbolic tangent activation functions a quadratic fit can be
obtained using a second-order Taylor series expansion of the activation
function $\mu(x)$ at the expansion point $x^{*}=0.5(l+u)$:
$\begin{split}&\mu(x)\approx\mu(x^{*})+\mu^{\prime}(x^{*})(x-x^{*})+0.5\,\mu^{\prime\prime}(x^{*})(x-x^{*})^{2}=\\\
&~{}\underbrace{0.5\,\mu^{\prime\prime}(x^{*})}_{a_{1}}x^{2}+\big{(}\underbrace{\mu^{\prime}(x^{*})-\mu^{\prime\prime}(x^{*})\,x^{*}}_{a_{2}}\big{)}x+\underbrace{\mu(x^{*})-\mu^{\prime}(x^{*})x^{*}+0.5\,\mu^{\prime\prime}(x^{*})\,{x^{*}}^{2}}_{a_{3}},\end{split}$
where the derivatives for sigmoid activations are
$\mu^{\prime}(x)=\sigma(x)(1-\sigma(x))$ and
$\mu^{\prime\prime}(x)=\sigma(x)(1-\sigma(x))(1-2\sigma(x))$, and the
derivatives for hyperbolic tangent activations are
$\mu^{\prime}(x)=1-\tanh(x)^{2}$ and
$\mu^{\prime\prime}(x)=-2\tanh(x)(1-\tanh(x)^{2})$. The Taylor series
expansion method is identical to the concept used in [21].

#### Linear Approximation

Since a linear function represents a special case of a quadratic function,
Alg. 1 can also be used in combination with linear approximations. Such
approximations are provided by the zonotope abstraction in [34]. Since closed-
form expressions for the bounds $\underline{d}$ and $\overline{d}$ of the
approximation error are already specified in [34], we can omit the error bound
computation described in Sec. 3.2 in this case. For ReLU activations we obtain
according to [34, Theorem 3.1]
$a_{1}=0,~{}~{}a_{2}=\frac{u}{u-l},~{}~{}a_{3}=\frac{-u\,l}{2(u-l)},~{}~{}\underline{d}=\frac{-u\,l}{2(u-l)},~{}~{}\overline{d}=\frac{u\,l}{2(u-l)},$
which results in the zonotope enclosure visualized in Fig. 1. For sigmoid and
hyperbolic tangent activations we obtain according to [34, Theorem 3.2]
$\begin{split}&a_{1}=0,~{}~{}a_{2}=\min(\mu^{\prime}(l),\mu^{\prime}(u)),~{}~{}a_{3}=0.5(\mu(u)+\mu(l)-a_{2}(u+l)),\\\
&\underline{d}=0.5(\mu(u)-\mu(l)-a_{2}(u-l)),~{}~{}\overline{d}=-0.5(\mu(u)-\mu(l)-a_{2}(u-l)),\end{split}$
where the derivatives of the sigmoid function and the hyperbolic tangent are
specified in the paragraph above.  Figure 3: Exemplary neural network with
ReLU activations (left) and the corresponding image enclosure computed with
polynomial zonotopes (right), where we use the approximation
$g(x)=0.25\,x^{2}+0.5\,x+0.25$ for the red neuron and the approximation
$g(x)=x$ for the blue neuron. We observed from experiments that for ReLU
activations the closed-form expression usually results in a tighter enclosure
of the image than polynomial regression. For sigmoid and hyperbolic tangent
activations, on the other hand, polynomial regression usually performs better
than the Taylor series expansion. It is also possible to combine multiple of
the methods described above by executing them in parallel and selecting the
one that results in the smallest approximation error
$[\underline{d},\overline{d}]$. Since the linear approximation does not
increase the number of generators, it represents an alternative to order
reduction when dealing with deep neural networks. Here, the development of a
method to decide automatically for which layers to use a linear and for which
a quadratic approximation is a promising direction for future research.

###  3.2 Bounding the Approximation Error

To obtain a sound enclosure we need to compute the difference between the
activation function $\mu(x)$ and the quadratic approximation
$g(x)=a_{1}\,x^{2}+a_{2}\,x+a_{3}$ from Sec. 3.1 on the interval $[l,u]$. In
particular, this corresponds to determining

$\underline{d}=\min_{x\in[l,u]}\underbrace{\mu(x)-a_{1}\,x^{2}-a_{2}\,x-a_{3}}_{d(x)}~{}~{}\text{and}~{}~{}\overline{d}=\max_{x\in[l,u]}\underbrace{\mu(x)-a_{1}\,x^{2}-a_{2}\,x-a_{3}}_{d(x)}.$

Depending on the type of activation function, we use different methods for
this.

#### Rectified Linear Unit (ReLU)

For ReLU activation functions we split the interval $[l,u]$ into the two
intervals $[l,0]$ and $[0,u]$ on which the activation function is constant and
linear, respectively. On the interval $[l,0]$ we have
$d(x)=-a_{1}\,x^{2}-a_{2}\,x-a_{3}$, and on the interval $[0,u]$ we have
$d(x)=-a_{1}\,x^{2}+(1-a_{2})\,x-a_{3}$. In both cases $d(x)$ is a quadratic
function whose maximum and minimum values are either located on the interval
boundary or at the point $x^{*}$ where the derivative of $d(x)$ is equal to
zero. The lower bound on $[l,0]$ is therefore given as
$\underline{d}=\min(d(l),d(x^{*}),d(0))$ if $x^{*}\in[l,0]$ and
$\underline{d}=\min(d(l),d(0))$ if $x^{*}\not\in[l,0]$, where
$x^{*}=-0.5\,a_{2}/a_{1}$. The upper bound as well as the bounds for $[0,u]$
are computed in a similar manner. Finally, the overall bounds are obtained by
taking the minimum and maximum of the bounds for the intervals $[l,0]$ and
$[0,u]$.

#### Sigmoid and Hyperbolic Tangent

Here our high-level idea is to sample the function $d(x)$ at points $x_{i}$
with distance $\Delta x$ distributed uniformly over the interval $[l,u]$. From
rough bounds for the derivative $d^{\prime}(x)$ we can then deduce how much
the function value between two sample points changes at most, which yields
tight bounds $\overline{d}_{b}\geq\overline{d}$ and
$\underline{d}_{b}\leq\underline{d}$. In particular, we want to choose the
sampling rate $\Delta x$ such that the bounds
$\overline{d}_{b},\underline{d}_{b}$ comply to a user-defined precision
$\delta>0$:
$\overline{d}+\delta\geq\overline{d}_{b}\geq\overline{d}~{}~{}\text{and}~{}~{}\underline{d}-\delta\leq\underline{d}_{b}\leq\underline{d}.$
(2) We observe that for both, sigmoid and hyperbolic tangent, the derivative
is globally bounded by $\mu^{\prime}(x)\in[0,\overline{\mu}]$, where
$\overline{\mu}=0.25$ for the sigmoid and $\overline{\mu}=1$ for the
hyperbolic tangent. In addition, it holds that the derivative of the quadratic
approximation $g(x)=a_{1}\,x^{2}+a_{2}\,x+a_{3}$ is bounded by
$g^{\prime}(x)\in[\underline{g},\overline{g}]$ on the interval $[l,u]$, where
$\underline{g}=\min(2a_{1}l+a_{2},2a_{1}u+a_{2})$ and
$\overline{g}=\max(2a_{1}l+a_{2},2a_{1}u+a_{2})$. As a consequence, the
derivative of the difference $d(x)=\mu(x)-g(x)$ is bounded by
$d^{\prime}(x)\in[-\overline{g},\overline{\mu}-\underline{g}]$. The value of
$d(x)$ can therefore at most change by $\pm\Delta d$ between two samples
$x_{i}$ and $x_{i+1}$, where $\Delta d=\Delta
x\max(|-\overline{g}|,|\overline{\mu}-\overline{g}|)$. To satisfy (2) we
require $\Delta d\leq\delta$, so that we have to choose the sampling rate as
$\Delta x\leq\delta/\max(|-\overline{g}|,|\overline{\mu}-\overline{g}|)$.
Finally, the bounds are computed by taking the maximum and minimum of all
samples: $\overline{d}_{b}=\max_{i}d(x_{i})+\delta$ and
$\underline{d}_{b}=\min_{i}d(x_{i})-\delta$. For our experiments we use a
precision of $\delta=0.001$.

##  4 Neural Network Controlled Systems

Reachable sets for neural network controlled systems can be computed
efficiently by combining our novel image enclosure approach for neural
networks with a reachability algorithm for nonlinear systems. We consider
general nonlinear systems

$\dot{x}(t)=f\big{(}x(t),u_{c}(x(t),t),w(t)\big{)},$ (3)

where $x\in\mathbb{R}^{n}$ is the system state,
$u_{c}:~{}\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{m}$ is a control law,
$w(t)\in\mathcal{W}\subset\mathbb{R}^{r}$ is a vector of uncertain
disturbances, and
$f:~{}\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{r}\to\mathbb{R}^{n}$
is a Lipschitz continuous function. For neural network controlled systems the
control law $u_{c}(x(t),t)$ is given by a neural network. Since neural network
controllers are usually realized as digital controllers, we consider the
sampled-data case where the control input is only updated at discrete times
$t_{0},t_{0}+\Delta t,t_{0}+2\Delta t,\dots,t_{F}$ and kept constant in
between. Here, $t_{0}$ is the initial time, $t_{F}$ is the final time, and
$\Delta t$ is the sampling rate. Without loss of generality, we assume from
now on that $t_{0}=0$ and $t_{F}$ is a multiple of $\Delta t$. The reachable
set is defined as follows:

###### Definition 3

(Reachable set) Let $\xi(t,x_{0},u_{c}(\cdot),w(\cdot))$ denote the solution
to (3) for initial state $x_{0}=x(0)$, control law $u_{c}(\cdot)$, and the
disturbance trajectory $w(\cdot)$. The reachable set for an initial set
$\mathcal{X}_{0}\subset\mathbb{R}^{n}$ and a disturbance set
$\mathcal{W}\subset\mathbb{R}^{r}$ is

$\mathcal{R}(t):=\big{\\{}\xi(t,x_{0},u_{c}(\cdot),w(\cdot))~{}\big{|}~{}x_{0}\in\mathcal{X}_{0},\forall
t^{*}\in[0,t]:~{}w(t^{*})\in\mathcal{W}\big{\\}}.$

Since the exact reachable set cannot be computed for general nonlinear
systems, we compute a tight enclosure instead. We exploit that the control
input is piecewise constant, so that the reachable set for each control cycle
can be computed using the extended system

$\begin{bmatrix}\dot{x}(t)\\\
\dot{u}(t)\end{bmatrix}=\begin{bmatrix}f(x(t),u(t),w(t))\\\
\mathbf{0}\end{bmatrix}$ (4)

together with the initial set $\mathcal{X}_{0}\times\mathcal{Y}$, where
$\mathcal{Y}$ is the image of $\mathcal{X}_{0}$ through the neural network
controller. The overall algorithm is specified in Alg. 2. Its high-level
concept is to loop over all control cycles, where in each cycle we first
compute the image of the current reachable set through the neural network
controller in Line 3. Next, the image is combined with the reachable set using
the Cartesian product in Line 4. This yields the initial set for the extended
system in (4), for which we compute the reachable set
$\widehat{\mathcal{R}}(t_{i+1})$ at time $t_{i+1}$ as well as the reachable
set $\widehat{\mathcal{R}}(\tau_{i})$ for the time interval $\tau_{i}$ in Line
6. While it is possible to use arbitrary reachability algorithms for nonlinear
systems, we apply the conservative polynomialization algorithm [2] since it
performs especially well in combination with polynomial zonotopes. Finally, in
Line 7, we project the reachable set back to the original system dimensions.

Algorithm 2 Reachable set for a neural network controlled system

Require: Nonlinear system $\dot{x}(t)=f(x(t),u_{c}(x(t),t),w(t))$, neural
network controller $u_{c}(x(t),t)$, initial set $\mathcal{X}_{0}$, disturbance
set $\mathcal{W}$, final time $t_{F}$, sampling rate $\Delta t$.

Ensure: Tight enclosure $\mathcal{R}\supseteq\mathcal{R}([0,t_{F}])$ of the
reachable set $\mathcal{R}([0,t_{F}])$.

1:$t_{0}\leftarrow 0$, $\mathcal{R}(t_{0})\leftarrow\mathcal{X}_{0}$

2:for $i\leftarrow 0$ to $t_{F}/\Delta t-1$ do (loop over all control cycles)

3: $\mathcal{Y}\leftarrow$ image of $\mathcal{R}(t_{i})$ through the neural
network controller using Alg. 1

4: $\widehat{\mathcal{R}}(t_{i})\leftarrow\mathcal{R}(t_{i})\times\mathcal{Y}$
(combine reachable set and input set using (7))

5: $t_{i+1}\leftarrow t_{i}+\Delta t$, $\tau_{i}\leftarrow[t_{i},t_{i+1}]$
(update time)

6: $\widehat{\mathcal{R}}(t_{i+1}),\widehat{\mathcal{R}}(\tau_{i})\leftarrow$
reachable set for extended system in (4) starting from
$\widehat{\mathcal{R}}(t_{i})$

7:
$\mathcal{R}(t_{i+1})\leftarrow[I_{n}~{}\mathbf{0}]\,\widehat{\mathcal{R}}(t_{i+1})$,
$\mathcal{R}(\tau_{i})\leftarrow[I_{n}~{}\mathbf{0}]\,\widehat{\mathcal{R}}(\tau_{i})$
(projection using (5))

8:end for

9:$\mathcal{R}\leftarrow\bigcup_{i=0}^{t_{F}/\Delta t-1}\mathcal{R}(\tau_{i})$
(reachable set for the whole time horizon)

##  5 Operations on Polynomial Zonotopes

Alg. 1 and Alg. 2 both require some special operations on polynomial
zonotopes, the implementation of which we present now. Given a polynomial
zonotope $\mathcal{PZ}=\langle c,G,G_{I},E\rangle_{PZ}\subset\mathbb{R}^{n}$,
a matrix $A\in\mathbb{R}^{o\times n}$, a vector $b\in\mathbb{R}^{o}$, and an
interval $\mathcal{I}=[l,u]\subset\mathbb{R}^{n}$, the affine map and the
Minkowski sum with an interval are given as

$\displaystyle A\,\mathcal{PZ}\oplus b=\langle Ac+b,AG,AG_{I},E\rangle_{PZ}$
(5) $\displaystyle\mathcal{PZ}\oplus\mathcal{I}=\langle
c+0.5(u+l),G,[G_{I}~{}0.5\,\text{diag}(u-l)],E\rangle_{PZ},$ (6)

which follows directly from [25, Prop. 8], [25, Prop. 9], and [1, Prop. 2.1].
For the Cartesian product used in Line 4 of Alg. 2 we can exploit the special
structure of the sets to calculate the Cartesian product of two polynomial
zonotopes $\mathcal{PZ}_{1}=\langle
c_{1},G_{1},G_{I,1},E_{1}\rangle_{PZ}\subset\mathbb{R}^{n}$ and
$\mathcal{PZ}_{2}=\langle
c_{2},[G_{2}~{}\widehat{G}_{2}],[G_{I,2}~{}\widehat{G}_{I,2}],[E_{1}~{}E_{2}]\rangle_{PZ}\subset\mathbb{R}^{o}$
as

$\mathcal{PZ}_{1}\times\mathcal{PZ}_{2}=\bigg{\langle}\begin{bmatrix}c_{1}\\\
c_{2}\end{bmatrix},\begin{bmatrix}G_{1}&\mathbf{0}\\\
G_{2}&\widehat{G}_{2}\end{bmatrix},\begin{bmatrix}G_{I,1}&\mathbf{0}\\\
G_{I,2}&\widehat{G}_{I,2}\end{bmatrix},[E_{1}~{}E_{2}]\bigg{\rangle}_{PZ}.$
(7)

In contrast to [25, Prop. 11], this implementation of the Cartesian product
explicitly preserves dependencies between the two sets, which is possible
since both polynomial zonotopes have identical dependent factors. Computing
the exact bounds of a polynomial zonotope in Line 6 of Alg. 1 would be
computationally infeasible, especially since this has to be done for each
neuron in the network. We therefore compute a tight enclosure of the bounds
instead, which can be done very efficiently:

###### Proposition 1

(Interval enclosure) Given a polynomial zonotope $\mathcal{PZ}=\langle
c,G,\linebreak[3]G_{I},E\rangle_{PZ}\subset\mathbb{R}^{n}$, an enclosing
interval can be computed as

$\mathcal{I}=[c+g_{1}-g_{2}-g_{3}-g_{4},c+g_{1}+g_{2}+g_{3}+g_{4}]\supseteq\mathcal{PZ}$

with

$\begin{split}&g_{1}=0.5\sum_{i\in\mathcal{H}}G_{(\cdot,i)},~{}g_{2}=0.5\sum_{i\in\mathcal{H}}|G_{(\cdot,i)}|,~{}g_{3}=\sum_{i\in\mathcal{K}}|G_{(\cdot,i)}|,~{}g_{4}=\sum_{i=1}^{q}|G_{I(\cdot,i)}|\\\
&~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{H}=\bigg{\\{}i~{}\bigg{|}~{}\prod_{j=1}^{p}\big{(}1-E_{(j,i)}\,\mathrm{mod}~{}2)\big{)}=1\bigg{\\}},~{}~{}\mathcal{K}=\\{1,\dots,h\\}\setminus\mathcal{H},\end{split}$

where $x\,\mathrm{mod}\,y$, $x,y\in\mathbb{N}_{0}$ is the modulo operation and
$\setminus$ denotes the set difference.

###### Proof

We first enclose the polynomial zonotope by a zonotope
$\mathcal{Z}\supseteq\mathcal{PZ}$ according to [25, Prop. 5], and then
compute an interval enclosure $\mathcal{I}\supseteq\mathcal{Z}$ of this
zonotope according to [1, Prop. 2.2].

The core operation for Alg. 1 is the computation of the image through a
quadratic function. While it is possible to obtain the exact image by
introducing new dependent factors, we compute a tight enclosure for
computational reasons:

###### Proposition 2

(Image quadratic function) Given a polynomial zonotope $\mathcal{PZ}=\langle
c,G,\linebreak[3]G_{I},E\rangle_{PZ}\subset\mathbb{R}$ and a quadratic
function $g(x)=a_{1}\,x^{2}+a_{2}\,x+a_{3}$ with
$a_{1},a_{2},a_{3},x\in\mathbb{R}$, the image of $\mathcal{PZ}$ through $g(x)$
can be tightly enclosed by

$\big{\\{}g(x)~{}\big{|}~{}x\in\mathcal{PZ}\big{\\}}\subseteq\langle
c_{q},G_{q},G_{I,q},E_{q}\rangle_{PZ}$

with

$\begin{split}&c_{q}=a_{1}c^{2}+a_{2}c+a_{3}+0.5\,a_{1}\sum_{i=1}^{q}G_{I(\cdot,i)}^{2},~{}~{}G_{q}=\big{[}(2a_{1}c+a_{2})G~{}~{}a_{1}\widehat{G}\big{]},\\\
&E_{q}=\big{[}E~{}~{}\widehat{E}\big{]},~{}~{}G_{I,q}=\big{[}(2a_{1}c+a_{2})G_{I}~{}~{}2a_{1}\overline{G}~{}~{}a_{1}\widecheck{G}\big{]},\end{split}$
(8)

where

$\begin{split}&\widehat{G}=\big{[}G^{2}~{}~{}2\,\widehat{G}_{1}~{}~{}\dots~{}~{}2\,\widehat{G}_{h-1}\big{]},~{}~{}\widehat{E}=\big{[}2\,E~{}~{}\widehat{E}_{1}~{}~{}\dots~{}~{}\widehat{E}_{h-1}\big{]},\\\
&\widehat{G}_{i}=\big{[}G_{(i)}G_{(i+1)}~{}~{}\dots~{}~{}G_{(i)}G_{(h)}\big{]},~{}~{}i=1,\dots,h-1,\\\
&\widehat{E}_{i}=\big{[}E_{(\cdot,i)}+E_{(\cdot,i+1)}~{}~{}\dots~{}~{}E_{(\cdot,i)}+E_{(\cdot,h)}\big{]},~{}~{}i=1,\dots,h-1,\\\
&\overline{G}=\big{[}G_{(1)}G_{I}~{}~{}\dots~{}~{}G_{(h)}G_{I}\big{]},~{}~{}\widecheck{G}=\big{[}0.5\,G_{I}^{2}~{}~{}2\,\widecheck{G}_{1}~{}~{}\dots~{}~{}2\,\widecheck{G}_{q-1}\big{]},\\\
&\widecheck{G}_{i}=\big{[}G_{I(i)}G_{I(i+1)}~{}~{}\dots~{}~{}G_{I(i)}G_{I(q)}\big{]},~{}~{}i=1,\dots,q-1,\end{split}$
(9)

and the squares in $G^{2}$ as well as $G_{I}^{2}$ are interpreted elementwise.

###### Proof

The proof is provided in Appendix A.

[width=0.98]./figures/image

Figure 4: Image enclosures computed with zonotopes (red), star sets (green),
Taylor models (purple), and polynomial zonotopes (blue) for randomly generated
neural networks with ReLU activations (left), sigmoid activations (middle),
and hyperbolic tangent activations (right). The exact image is shown in gray.

##  6 Numerical Examples

We now demonstrate the performance of our approach for image computation,
open-loop neural network verification, and reachability analysis of neural
network controlled systems. If not stated otherwise, computations are carried
out in MATLAB on a 2.9GHz quad-core i7 processor with 32GB memory. We
integrated our implementation into CORA [3] and published a repeatability
package-1-1-1https://codeocean.com/capsule/8237552/tree/v1.

#### Image Enclosure

First, we demonstrate how our approach captures the non-convexity of the image through a neural network. For visualization purposes we use the deliberately simple example of randomly generated neural networks with two inputs, two outputs, and one hidden layer consisting of 50 neurons. The initial set is $\mathcal{X}_{0}=[-1,1]\times[-1,1]$. We compare our polynomial-zonotope-based approach with the zonotope abstraction in [34], the star set approach in [40] using the triangle relaxation, and the Taylor model abstraction in [21]. While our approach and the zonotope abstraction are applicable to all types of activation functions, the star set approach is restricted to ReLU activations and the Taylor model abstraction is limited to sigmoid and hyperbolic tangent activations. The resulting image enclosures are visualized in Fig. 4. While using zonotopes or star sets only yields a convex over-approximation, polynomial zonotopes are able to capture the non-convexity of the image and therefore provide a tighter enclosure. While Taylor models also capture the non-convexity of the image to some extent they are less precise than polynomial zonotopes, which can be explained as follows: 1) The zonotopic remainder of polynomial zonotopes prevents the rapid remainder growth observed for Taylor models, and 2) the quadratic approximation obtained with polynomial regression used for polynomial zonotopes is usually more precise than the Taylor series expansion used for Taylor models.  Table 1: Computation times00footnotemark: 0 in seconds for different verification tools on a small but representative excerpt of network-specification combinations of the ACAS Xu benchmark. The symbol - indicates that the tool failed to verify the specification. Net. | Spec. |  Cgdtest |  CROWN |  Debona |  ERAN |  Marabou |  MN-BaB |  nnenum |  nnv |  NV.jl |  oval |  RPM |  venus2 |  VeriNet |  Poly. zono.
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1.9 | 1 | 0.37 | 1.37 | 111 | 3.91 | 0.66 | 48.7 | 0.41 | - | 1.44 | 0.71 | - | 0.53 | 0.55 | 0.31
2.3 | 4 | - | 0.95 | 1.78 | 1.91 | 0.57 | 12.2 | 0.06 | - | - | 0.97 | - | 0.46 | 0.17 | 0.16
3.5 | 3 | 0.41 | 0.37 | 1.15 | 1.85 | 0.61 | 6.17 | 0.05 | - | - | 0.58 | 34.1 | 0.42 | 0.25 | 0.32
4.5 | 4 | - | 0.35 | 0.20 | 1.82 | 0.61 | 5.57 | 0.08 | 0.24 | - | 0.48 | - | 0.42 | 0.21 | 0.16
5.6 | 3 |  0.38 |  0.63 |  2.27 |  1.82 |  0.66 |  6.51 |  0.08 |  - |  - |  0.52 |  40.6 |  0.48 |  0.37 |  0.43

#### Open-Loop Neural Network Verification

For open-loop neural network verification the task is to verify that the image of the neural network satisfies certain specifications that are typically given by linear inequality constraints. We examine the ACAS Xu benchmark from the 2021 and 2022 VNN competition [4, 29] originally proposed in [22, Sec. 5], which features neural networks that provide turn advisories for an aircraft to avoid collisions. All networks consist of 6 hidden layers with 50 ReLU neurons per layer. For a fair comparison we performed the evaluation on the same machine that was used for the VNN competition. To compute the image through the neural networks with polynomial zonotopes, we apply a quadratic approximation obtained by polynomial regression for the first two layers, and a linear approximation in the remaining layers. Moreover, we recursively split the initial set to obtain a complete verifier. The comparison with the other tools that participated in the VNN competition shown in Tab. 1 demonstrates that for some verification problems polynomial zonotopes are about as fast as the best tool in the competition. 00footnotetext: Times taken from https://github.com/stanleybak/vnncomp2021_results and https://github.com/ChristopherBrix/vnncomp2022_results. Table 2: Computation times111Computation times taken from [33, Tab. 1] for Sherlock and JuliaReach, from [21, Tab. 2] for Verisig, Verisig 2.0, and ReachNN*, and from [19, Tab. 1] for POLAR. in seconds for reachability analysis of neural network controlled systems considering different tools and benchmarks. The dimension, the number of hidden layers, and the number of neurons in each layer is specified in parenthesis for each benchmark, where $a=100$, $b=5$ for ReLU activation functions, and $a=20$, $b=3$ otherwise. The symbol - indicates that the tool failed to verify the specification. | ReLU | sigmoid | hyp. tangent
---|---|---|---
|  Sherlock |  JuliaReach |  Poly. zono. |  Verisig |  Verisig 2.0 |  ReachNN* |  POLAR |  Poly. zono. |  Verisig |  Verisig 2.0 |  ReachNN* |  POLAR |  Poly. zono.
B1 (2, 2, 20) |  |  |  | - | 49 | 69 | 23 | 2 | - | 48 | - | 25 | 8
B2 (2, 2, 20) |  |  |  | 12 | 8 | 32 | 10 | 1 | - | - | - | 3 | -
B3 (2, 2, 20) |  |  |  | 98 | 47 | 130 | 37 | 3 | 98 | 43 | 128 | 38 | 3
B4 (3, 2, 20) |  |  |  | 24 | 12 | 20 | 4 | 1 | 23 | 11 | 20 | 4 | 1
B5 (3, 3, 100) |  |  |  |  196 | 1063 | 31 | 25 | 2 | - | 168 | - | 31 | 2
TORA (4, 3, $a$) | 30 | 2040 | 13 | 136 | 83 | 13402 |  | 1 | 134 | 70 | 2524 |  | 1
ACC (6, $b$, 20) | 4 | 1 | 2 |  |  |  |  |  | - | 1512 | - | 312 | 2
Unicycle (3, 1, 500) | 526 | 93 | 3 |  |  |  |  |  |  |  |  |  |
Airplane (12, 3, 100) | - | 29 | 7 |  |  |  |  |  |  |  |  |  |
Sin. Pend. (2, 2, 25) | 1 | 1 | 1 |  |  |  |  |  |  |  |  |  |

#### Neural Network Controlled Systems

The main application of our approach is reachability analysis of neural
network controlled systems, for which we now compare the performance to other
state-of-the-art tools. For a fair comparison we focus on published results
for which the authors of the tools tuned the algorithm settings by themselves.
In particular, we examine the benchmarks from [33] featuring ReLU neural
network controllers, and the benchmarks from [21] containing sigmoid and
hyperbolic tangent neural network controllers. The goal for all benchmarks is
to verify that the system reaches a goal set or avoids an unsafe region. As
the computation times shown in Tab. 1 demonstrate, our polynomial-zonotope-
based approach is for all but two benchmarks significantly faster than all
other state-of-the-art tools, mainly since it avoids all major bottlenecks
observed for the other tools: The polynomial approximations of the overall
network used by Sherlock and ReachNN* are often imprecise, JuliaReach loses
dependencies when enclosing Taylor models by zonotopes, Verisig is quite slow
since the nonlinear system used to represent the neural network is high-
dimensional, and Verisig 2.0 and POLAR suffer from the rapid remainder growth
observed for Taylor models.

##  7 Conclusion

We introduced a novel approach for computing tight non-convex enclosures of
images through neural networks with ReLU, sigmoid, and hyperbolic tangent
activation functions. Since we represent sets with polynomial zonotopes, all
required calculations can be realized using simple matrix operations only,
which makes our algorithm very efficient. While our proposed approach can also
be applied to open-loop neural network verification, its main application is
reachability analysis of neural network controlled systems. There, polynomial
zonotopes enable the preservation of dependencies between the reachable set
and the set of control inputs, which results in very tight enclosures of the
reachable set. As we demonstrated on various numerical examples, our
polynomial-zonotope-based approach consequently outperforms all other state-
of-the-art methods for reachability analysis of neural network controlled
systems.

Acknowledgements. We gratefully acknowledge the financial support from the
project justITSELF funded by the European Research Council (ERC) under grant
agreement No 817629, from DIREC - Digital Research Centre Denmark, and from
the Villum Investigator Grant S4OS. In addition, this material is based upon
work supported by the Air Force Office of Scientific Research and the Office
of Naval Research under award numbers FA9550-19-1-0288, FA9550-21-1-0121,
FA9550-23-1-0066 and N00014-22-1-2156. Any opinions, findings, and conclusions
or recommendations expressed in this material are those of the authors and do
not necessarily reflect the views of the United States Air Force or the United
States Navy.

## References

  * [1] Althoff, M.: Reachability Analysis and its Application to the Safety Assessment of Autonomous Cars. Ph.D. thesis, Technical University of Munich (2010)
  * [2] Althoff, M.: Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets. In: Proc. of the Int. Conf. on Hybrid Systems: Computation and Control. pp. 173–182 (2013)
  * [3] Althoff, M.: An introduction to CORA 2015. In: Proc. of the Int. Workshop on Applied Verification for Continuous and Hybrid Systems. pp. 120–151 (2015)
  * [4] Bak, S., Liu, C., Johnson, T.: The second international verification of neural networks competition (VNN-COMP 2021): Summary and results. arXiv:2109.00498 (2021)
  * [5] Beaufays, F., Abdel-Magid, Y., Widrow, B.: Application of neural networks to load-frequency control in power systems. Neural Networks 7(1), 183–194 (1994)
  * [6] Bogomolov, S., et al.: JuliaReach: A toolbox for set-based reachability. In: Proc. of the Int. Conf. on Hybrid Systems: Computation and Control. pp. 39–44 (2019)
  * [7] Bunel, R., et al.: Branch and bound for piecewise linear neural network verification. Journal of Machine Learning Research 21(42) (2020)
  * [8] Cheng, C.H., Nührenberg, G., Ruess, H.: Maximum resilience of artificial neural networks. In: Proc. of the Int. Symposium on Automated Technology for Verification and Analysis. pp. 251–268 (2017)
  * [9] Christakou, C., Vrettos, S., Stafylopatis, A.: A hybrid movie recommender system based on neural networks. International Journal on Artificial Intelligence Tools 16(5), 771–792 (2007)
  * [10] Clavière, A., et al.: Safety verification of neural network controlled systems. In: Proc. of the Int. Conf. on Dependable Systems and Networks. pp. 47–54 (2021)
  * [11] David, O.E., Netanyahu, N.S., Wolf, L.: DeepChess: End-to-end deep neural network for automatic learning in chess. In: Proc. of the Int. Conf. on Artificial Neural Networks. pp. 88–96 (2016)
  * [12] Dutta, S., et al.: Learning and verification of feedback control systems using feedforward neural networks. In: Proc. of the Int. Conf. on Analysis and Design of Hybrid Systems. pp. 151–156 (2018)
  * [13] Dutta, S., et al.: Sherlock–A tool for verification of neural network feedback systems. In: Proc. of the Int. Conf. on Hybrid Systems: Computation and Control. pp. 262–263 (2019)
  * [14] Dutta, S., Chen, X., Sankaranarayanan, S.: Reachability analysis for neural feedback systems using regressive polynomial rule inference. In: Proc. of the Int. Conf. on Hybrid Systems: Computation and Control. pp. 157–168 (2019)
  * [15] Ehlers, R.: Formal verification of piece-wise linear feed-forward neural networks. In: Proc. of the Int. Symposium on Automated Technology for Verification and Analysis. pp. 269–286 (2017)
  * [16] Fan, J., Huang, et al.: ReachNN*: A tool for reachability analysis of neural-network controlled systems. In: Proc. of the Int. Symposium on Automated Technology for Verification and Analysis. pp. 537–542 (2020)
  * [17] Goubault, E., Putot, S.: RINO: Robust inner and outer approximated reachability of neural networks controlled systems. In: Proc. of the Int. Conf. on Computer Aided Verification. pp. 511–523 (2022)
  * [18] Huang, C., et al.: ReachNN: Reachability analysis of neural-network controlled systems. Transactions on Embedded Computing Systems 18(5s) (2019)
  * [19] Huang, C., et al.: POLAR: A polynomial arithmetic framework for verifying neural-network controlled systems. In: Proc. of the Int. Symposium on Automated Technology for Verification and Analysis. pp. 414–430 (2022)
  * [20] Ivanov, R., et al.: Verisig: Verifying safety properties of hybrid systems with neural network controllers. In: Proc. of the Int. Conf. on Hybrid Systems: Computation and Control. pp. 169–178 (2019)
  * [21] Ivanov, R., et al.: Verisig 2.0: Verification of neural network controllers using Taylor model preconditioning. In: Proc. of the Int. Conf. on Computer Aided Verification. pp. 249–262 (2021)
  * [22] Katz, G., et al.: Reluplex: An efficient SMT solver for verifying deep neural networks. In: Proc. of the Int. Conf. on Computer Aided Verification. pp. 97–117 (2017)
  * [23] Khan, S., et al.: Facial recognition using convolutional neural networks and implementation on smart glasses. In: Proc. of the Int. Conf. on Information Science and Communication Technology (2019), Article 19
  * [24] Khedr, H., Ferlez, J., Shoukry, Y.: PEREGRiNN: Penalized-relaxation greedy neural network verifier. In: Proc. of the Int. Conf. on Computer Aided Verification. pp. 287–300 (2021)
  * [25] Kochdumper, N., Althoff, M.: Sparse polynomial zonotopes: A novel set representation for reachability analysis. Transactions on Automatic Control 66(9), 4043–4058 (2021)
  * [26] Makino, K., Berz, M.: Taylor models and other validated functional inclusion methods. Int. Journal of Pure and Applied Mathematics 4(4), 379–456 (2003)
  * [27] Mukherjee, D., et al.: A survey of robot learning strategies for human-robot collaboration in industrial settings. Robotics and Computer-Integrated Manufacturing 73 (2022)
  * [28] Müller, M.N., et al.: PRIMA: Precise and general neural network certification via multi-neuron convex relaxations. Proceedings on Programming Languages 1 (2022), Article 43
  * [29] Müller, M.N., et al.: The third international verification of neural networks competition (VNN-COMP 2022): Summary and results. arXiv preprint arXiv:2212.10376 (2022)
  * [30] Pulina, L., Tacchella, A.: Challenging SMT solvers to verify neural networks. AI Communications 25(2), 117–135 (2012)
  * [31] Raghunathan, A., Steinhardt, J., Liang, P.: Semidefinite relaxations for certifying robustness to adversarial examples. In: Proc. of the Int. Conf. on Neural Information Processing Systems. pp. 10900–10910 (2018)
  * [32] Riedmiller, M., Montemerlo, M., Dahlkamp, H.: Learning to drive a real car in 20 minutes. In: Proc. of the Int. Conf. on Frontiers in the Convergence of Bioscience and Information Technologies. pp. 645–650 (2007)
  * [33] Schilling, C., Forets, M., Guadalupe, S.: Verification of neural-network control systems by integrating Taylor models and zonotopes. In: Proc. of the AAAI Conf. on Artificial Intelligence. pp. 8169–8177 (2022)
  * [34] Singh, G., et al.: Fast and effective robustness certification. In: Proc. of the Int. Conf. on Advances in Neural Information Processing Systems (2018)
  * [35] Singh, G., et al.: An abstract domain for certifying neural networks. Proceedings on Programming Languages 3 (2019), Article 41
  * [36] Singh, G., et al.: Beyond the single neuron convex barrier for neural network certification. In: Proc. of the Int. Conf. on Advances in Neural Information Processing Systems (2019)
  * [37] Tjeng, V., Xiao, K.Y., Tedrake, R.: Evaluating robustness of neural networks with mixed integer programming. In: Proc. of the Int. Conf. on Learning Representations (2019)
  * [38] Tran, H.D., et al.: Parallelizable reachability analysis algorithms for feed-forward neural networks. In: Proc. of the Int. Conf. on Formal Methods in Software Engineering. pp. 51–60 (2019)
  * [39] Tran, H.D., et al.: Safety verification of cyber-physical systems with reinforcement learning control. Transactions on Embedded Computing Systems 18(5s) (2019), Article 105
  * [40] Tran, H.D., et al.: Star-based reachability analysis of deep neural networks. In: Proc. of the Int. Symposium on Formal Methods. pp. 670–686 (2019)
  * [41] Tran, H.D., et al.: NNV: The neural network verification tool for deep neural networks and learning-enabled cyber-physical systems. In: Proc. of the Int. Conf. on Computer Aided Verification. pp. 3–17 (2020)
  * [42] Vincent, J.A., Schwager, M.: Reachable Polyhedral Marching (RPM): A safety verification algorithm for robotic systems with deep neural network components. In: Proc. of the Int. Conf. on Robotics and Automation. pp. 9029–9035 (2021)
  * [43] Wang, S., et al.: Formal security analysis of neural networks using symbolic intervals. In: Proc. of the USENIX Security Symposium. pp. 1599–1614 (2018)
  * [44] Wang, S., et al.: Beta-CROWN: Efficient bound propagation with per-neuron split constraints for neural network robustness verification. In: Proc. of the Int. Conf. on Neural Information Processing Systems (2021)
  * [45] Weng, L., et al.: Towards fast computation of certified robustness for ReLU networks. In: Proc. of the Int. Conf. on Machine Learning. pp. 5276–5285 (2018)
  * [46] Xiang, W., et al.: Reachable set estimation and safety verification for piecewise linear systems with neural network controllers. In: Proc. of the American Control Conf. pp. 1574–1579 (2018)
  * [47] Xu, K., et al.: Fast and complete: Enabling complete neural network verification with rapid and massively parallel incomplete verifiers. In: Proc. of the Int. Conf. on Learning Representations (2021)
  * [48] Yang, X., et al.: Reachability analysis of deep ReLU neural networks using facet-vertex incidence. In: Proc. of the Int. Conference on Hybrid Systems: Computation and Control (2021), Article 18
  * [49] Zhang, H., et al.: Efficient neural network robustness certification with general activation functions. In: Proc. of the Int. Conf. on Neural Information Processing Systems. pp. 4944–4953 (2018)

## Appendix A

We now provide the proof for Prop. 2. According to Def. 2, the one-dimensional
polynomial zonotope $\mathcal{PZ}=\langle c,G,G_{I},E\rangle_{PZ}$ is defined
as

$\begin{split}\mathcal{PZ}&=\bigg{\\{}c+\underbrace{\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\bigg{)}G_{(i)}}_{d(\alpha)}+\underbrace{\sum_{j=1}^{q}\beta_{j}G_{I(j)}}_{z(\beta)}~{}\bigg{|}~{}\alpha_{k},\beta_{j}\in[-1,1]\bigg{\\}}\\\
&=\big{\\{}c+d(\alpha)+z(\beta)~{}\big{|}~{}\alpha,\beta\in[-\mathbf{1},\mathbf{1}]\big{\\}},\end{split}$
(10)

where $\alpha=[\alpha_{1}~{}\dots~{}\alpha_{p}]^{T}$ and
$\beta=[\beta_{1}~{}\dots~{}\beta_{q}]^{T}$. To compute the image through the
quadratic function $g(x)$ we require the expressions $d(\alpha)^{2}$,
$d(\alpha)z(\beta)$, and $z(\beta)^{2}$, which we derive first. For
$d(\alpha)^{2}$ we obtain

$\begin{split}d(\alpha)^{2}&=\bigg{(}\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\bigg{)}G_{(i)}\bigg{)}\bigg{(}\sum_{j=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,j)}}\bigg{)}G_{(j)}\bigg{)}\\\
&=\sum_{i=1}^{h}\sum_{j=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}+E_{(k,j)}}\bigg{)}G_{(i)}G_{(j)}\\\
&=\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{2E_{(k,i)}}\bigg{)}G_{(i)}^{2}+\sum_{i=1}^{h-1}\sum_{j=i+1}^{h}\bigg{(}\prod_{k=1}^{p}\underbrace{\alpha_{k}^{E_{(k,i)}+E_{(k,j)}}}_{\alpha_{k}^{\widehat{E}_{i(k,j)}}}\bigg{)}2\underbrace{G_{(i)}G_{(j)}}_{\widehat{G}_{i(j)}}\\\
&\overset{\begin{subarray}{c}\eqref{eq:quadMap2}\\\
\vspace{-2pt}\end{subarray}}{=}\sum_{i=1}^{h(h+1)/2}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{\widehat{E}_{(k,i)}}\bigg{)}\widehat{G}_{(i)},\end{split}$
(11)

for $d(\alpha)z(\beta)$ we obtain

$\begin{split}d(\alpha)z(\beta)&=\bigg{(}\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\bigg{)}G_{(i)}\bigg{)}\bigg{(}\sum_{j=1}^{q}\beta_{j}G_{I(j)}\bigg{)}\\\
&=\sum_{i=1}^{h}\sum_{j=1}^{q}\underbrace{\bigg{(}\beta_{j}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\bigg{)}}_{\beta_{q+(i-1)h+j}}G_{(i)}G_{I(j)}\overset{\begin{subarray}{c}\eqref{eq:quadMap2}\\\
\vspace{-2pt}\end{subarray}}{=}\sum_{i=1}^{hq}\beta_{q+i}\,\overline{G}_{(i)},\end{split}$
(12)

and for $z(\beta)^{2}$ we obtain

$\displaystyle z(\beta)^{2}$
$\displaystyle=\bigg{(}\sum_{i=1}^{q}\beta_{i}G_{I(i)}\bigg{)}\bigg{(}\sum_{j=1}^{q}\beta_{j}G_{I(j)}\bigg{)}=\sum_{i=1}^{q}\sum_{j=1}^{q}\beta_{i}\beta_{j}\,G_{I(i)}G_{I(j)}$
$\displaystyle=\sum_{i=1}^{q}\beta_{i}^{2}G_{I(i)}^{2}+\sum_{i=1}^{q-1}\sum_{j=i+1}^{q}\beta_{i}\beta_{j}\,2\,G_{I(i)}G_{I(j)}$
$\displaystyle=0.5\sum_{i=1}^{q}G_{I(i)}^{2}+\sum_{i=1}^{q}\underbrace{(2\beta_{i}^{2}-1)}_{\beta_{(h+1)q+i}}0.5\,G_{I(i)}^{2}+\sum_{i=1}^{q-1}\sum_{j=i+1}^{q}\underbrace{\beta_{i}\beta_{j}}_{\beta_{a(i,j)}}2\,\underbrace{G_{I(i)}G_{I(j)}}_{\widecheck{G}_{i(j)}}$
(13) $\displaystyle\overset{\begin{subarray}{c}\eqref{eq:quadMap2}\\\
\vspace{-2pt}\end{subarray}}{=}0.5\sum_{i=1}^{q}G_{I(i)}^{2}+\sum_{i=1}^{q(q+1)/2}\beta_{(h+1)q+i}\,\widecheck{G}_{(i)},$

where the function $a(i,j)$ maps indices $i,j$ to a new index:

$a(i,j)=(h+2)q+j-i+\sum_{k=1}^{i-1}q-k.$

In (12) and (13), we substituted the expressions
$\beta_{j}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}$, $2\beta_{i}^{2}-1$, and
$\beta_{i}\beta_{j}$ containing polynomial terms of the independent factors
$\beta$ by new independent factors, which results in an enclosure due to the
loss of dependency. The substitution is possible since

$\beta_{j}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\in[-1,1],~{}~{}2\beta_{i}^{2}-1\in[-1,1],~{}~{}\text{and}~{}~{}\beta_{i}\beta_{j}\in[-1,1].$

Finally, we obtain for the image

$\begin{split}&\big{\\{}g(x)~{}\big{|}~{}x\in\mathcal{PZ}\big{\\}}=\big{\\{}a_{1}\,x^{2}+a_{2}\,x+a_{3}~{}\big{|}~{}x\in\mathcal{PZ}\big{\\}}\overset{\begin{subarray}{c}\eqref{eq:defPolyZono}\\\
\vspace{-2pt}\end{subarray}}{=}\\\ &~{}\\\
&\big{\\{}a_{1}(c+d(\alpha)+z(\beta))^{2}+a_{2}(c+d(\alpha)+z(\beta))+a_{3}~{}\big{|}~{}\alpha,\beta\in[-\mathbf{1},\mathbf{1}]\big{\\}}=\\\
&~{}\\\
&\big{\\{}a_{1}c^{2}+a_{2}c+a_{3}+(2a_{1}c+a_{2})d(\alpha)+a_{1}d(\alpha)^{2}\\\
&~{}~{}+(2a_{1}c+a_{2})z(\beta)+2a_{1}d(\alpha)z(\beta)+a_{1}z(\beta)^{2}~{}\big{|}~{}\alpha,\beta\in[-\mathbf{1},\mathbf{1}]\big{\\}}\overset{\begin{subarray}{c}\eqref{eq:expression1},\eqref{eq:expression2},\eqref{eq:expression3}\\\
\vspace{-2pt}\end{subarray}}{\subseteq}\\\ &~{}\\\
&\bigg{\langle}a_{1}c^{2}+a_{2}c+a_{3}+0.5\,a_{1}\sum_{i=1}^{q}G_{I}^{2},\big{[}(2a_{1}c+a_{2})G~{}~{}a_{1}\widehat{G}\big{]},\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\big{[}(2a_{1}c+a_{2})G_{I}~{}~{}2a_{1}\overline{G}~{}~{}a_{1}\widecheck{G}\big{]},\big{[}E~{}~{}\widehat{E}\big{]}\bigg{\rangle}_{PZ}\overset{\begin{subarray}{c}\eqref{eq:quadMap1}\\\
\vspace{-2pt}\end{subarray}}{=}\langle
c_{q},G_{q},G_{I,q},E_{q}\rangle_{PZ},\end{split}$

which concludes the proof.


# Empowering Healthcare through Privacy-Preserving MRI Analysis

Al Amin1, Kamrul Hasan1, Saleh Zein-Sabatto1, Deo Chimba1, Liang Hong1, Imtiaz
Ahmed2, Tariqul Islam3 1 Tennessee State University, Nashville, TN, USA
2Howard University, Washington, DC, USA
3 Syracuse University, Syracuse, NY, USA
Email: $\\{$aamin2, mhasan1, mzein, dchimba<EMAIL_ADDRESS>$\\{$imtiaz.ahmed$\\}$ @howard.edu, $\\{$mtislam$\\}$ @syr.edu

###### Abstract

In the healthcare domain, Magnetic Resonance Imaging (MRI) assumes a pivotal
role, as it employs Artificial Intelligence (AI) and Machine Learning (ML)
methodologies to extract invaluable insights from imaging data. Nonetheless,
the imperative need for patient privacy poses significant challenges when
collecting data from diverse healthcare sources. Consequently, the Deep
Learning (DL) communities occasionally face difficulties detecting rare
features. In this research endeavor, we introduce the Ensemble-Based Federated
Learning (EBFL) Framework, an innovative solution tailored to address this
challenge. The EBFL framework deviates from the conventional approach by
emphasizing model features over sharing sensitive patient data. This unique
methodology fosters a collaborative and privacy-conscious environment for
healthcare institutions, empowering them to harness the capabilities of a
centralized server for model refinement while upholding the utmost data
privacy standards. Conversely, a robust ensemble architecture boasts potent
feature extraction capabilities, distinguishing itself from a single DL model.
This quality makes it remarkably dependable for MRI analysis. By harnessing
our groundbreaking EBFL methodology, we have achieved remarkable precision in
the classification of brain tumors, including glioma, meningioma, pituitary,
and non-tumor instances, attaining a precision rate of 94% for the Global
model and an impressive 96% for the Ensemble model. Our models underwent
rigorous evaluation using conventional performance metrics such as Accuracy,
Precision, Recall, and F1 Score. Integrating DL within the Federated Learning
(FL) framework has yielded a methodology that offers precise and dependable
diagnostics for detecting brain tumors.

###### Index Terms:

Federated Learning (FL), Maximum Voting Classifier (Ensemble), Data privacy,
Intelligent Healthcare System, Health.

##  I Introduction

Identifying a brain tumor is a critical and life-altering event for patients.
Globally, approximately 190,000 cases of both primary and secondary brain
tumors are diagnosed each year, impacting individuals of all age groups [1].
Physicians extensively utilize magnetic resonance (MR) imaging to ascertain
the presence and characteristics of tumors [2, 3, 4]. Therefore, researchers
diligently concentrate on AI, ML, and DL frameworks to extract meaningful
features from MR images [2, 5, 6, 7]. Nonetheless, achieving enhanced
performance with DL architecture necessitates a substantial volume of data [8,
9, 10, 11, 12]. The exchange of patient data among various healthcare
providers is complicated by the need to protect patient confidentiality [13,
14, 15]. This situation makes it difficult to gather more data from different
hospitals [16, 17]. In this research, our focus centers on the precise
identification of various brain tumor types, including glioma, meningioma,
pituitary tumors, and non-tumorous conditions, through the utilization of MR
images. To achieve this objective while respecting patient data privacy, we
introduce the novel EBFL framework. This innovative approach enables
individual hospitals to train their data with AI tools, sharing their feature-
extracted trained models at a central server aggregation point without needing
patient data exchange. The main server model leverages potent feature
extraction techniques, amalgamating insights from diverse hospitals.
Additionally, we employ an ensemble model, incorporating multiple DL
architectures, to counter model overfitting and harness various feature
extraction methods inherent to distinct DL architectures. This comprehensive
strategy enhances our diagnostic algorithm’s reliability and robustness,
surpassing a single algorithm’s capabilities.

  * •

First, the research introduces a robust Federated Learning (FL) model
incorporating a maximum voting classifier (Ensemble techniques) based on DL
architecture designed for a privacy-preserving hybrid paradigm. This
integration is intended to enhance diagnostic reliability and precision,
providing healthcare practitioners with improved tools for detecting brain
tumors efficiently.

  * •

Second, our research introduces the EBFL Framework. This novel approach
upholds patient data confidentiality while enabling multiple healthcare
institutions to benefit from a centralized diagnostic model. This framework
allows for synthesizing a global model’s intelligence, with participating
entities contributing to and deriving insights from the model without
exchanging sensitive data.

  * •

Finally, we developed a robust ensemble mode for training local data by
combining multiple Transfer Learning architectures MobileNetV2, VGG16, and
VGG19.

##  II Related Works

DL architectures are crucial in healthcare research, enhancing MRI analysis
for applications such as COVID-19 classification [18, 19, 20, 21], spinal cord
disorders [22, 23, 24], cancer detection and staging [25, 26, 27],
musculoskeletal disorders [28, 29], abdominal and pelvic disorders [30], and
brain and neurological disorders [31, 32, 33]. However, Brain tumor
recognition is a challenging task. A notable example of advancements in brain
tumor classification, specifically targeting glioma, meningioma, and pituitary
tumors, is the development of DeepTumorNet, a hybrid DL model. This model
innovatively adapts the GoogleNet architecture, a CNN variant. In creating
DeepTumorNet, the original architecture was modified by removing the last five
layers of GoogleNet and integrating fifteen new layers in their place. This
adaptation has led to remarkable outcomes. However, it’s important to note the
limitations of this approach. DeepTumorNet requires a substantial amount of
data for training and is characterized by its computational complexity, which
could pose challenges in practical healthcare applications due to patient
privacy[34]. Another research group applies a deep CNN for brain tumor
classification from MR images. A notable limitation is the small size of the
dataset used for training, which can lead to overfitting and reduced
generalizability of the model. Additionally, the high variance in MR images
due to different scanning protocols and machines poses a challenge for
accurate classification. The research suggests the need for larger,
standardized datasets and improved models for better diagnostic accuracy in
medical imaging [35]. A new set of researchers applied a two-channel deep
neural network (DNN) model for brain tumor classification, utilizing pre-
trained CNNs like InceptionResNetV2 and Xception for feature extraction. The
model emphasized tumor-affected regions using an attention mechanism and
achieved strong generalization on benchmark datasets. A limitation of this
research is its focus on only two specific pre-trained models. Future work
plans to explore additional pre-training models and develop a new method for
abstracting features from convolution blocks [36]. A diverse research team
developed a method for brain tumor classification that combines an ensemble of
deep features with machine learning classifiers. This approach utilizes
transfer learning with various pre-trained DCNN to extract features from brain
MR images, which are then evaluated by multiple classifiers, including a
notably effective support vector machine (SVM) with radial basis function
(RBF) kernel, particularly for large datasets [37]. Another consortium of
researchers used a methodology involving data acquisition, preprocessing,
feature extraction using texture features, feature selection, and
classification with a Support Vector Machine (SVM) ensemble model. The
limitation was the focus on intermediate-stage tumor data, which might impact
the model’s effectiveness in early-stage tumor classification. This is an area
for future improvement, along with enhancing the robustness of the SVM learner
and voting algorithm [38]. An alternative scientific group mentions various
applications of FL in medical imaging, like MRI reconstruction with
decentralized training of generative image priors, distinguishing healthy and
cancerous brain tissues, explicitly highlighting its use in brain tumor
segmentation and classification. The limitations of FL in this research area,
such as communication challenges, data bias, and the need for performance
improvement strategies, are also addressed [39]. In this research, we have
successfully addressed the limitations of existing FL applications in brain
tumor classification by developing an innovative ensemble model. This model
synergistically combines three DL architectures: MobileNetV2, VGG16, and
VGG19. Each architecture employs distinct feature extraction techniques, which
are particularly adept at handling small datasets with rare features. This
approach enhances performance and effectively mitigates the risk of model
overfitting. Furthermore, our implementation of federated learning plays a
crucial role in preserving patient privacy, as it focuses on feature analysis
without necessitating the sharing of raw data. Significantly, we have also
introduced a novel framework that developed communication among healthcare
industries.

##  III Methodology

The EBFL framework, as shown in Figure 1, integrates multiple DL algorithms
for MRI image analysis, using majority voting and feature sharing at the
aggregation point to enhance global model robustness and ensure data privacy,
effectively addressing the challenges of data scarcity.

Figure 1: Ensemble Based Federated Learning (EBFL) framework for analyzing MR
Images

###  III-A EBFL based MR images System Description

The proposed EBFL framework for brain tumor MRI analysis is structured into
three sequential phases, as illustrated in Figure 1. Initially, the framework
decentralizes the comprehensive brain tumor MRI dataset across four distinct
clients. Subsequently, each client independently trains an ensemble of models
on their respective server, leveraging the strength of combined algorithms for
robust feature extraction. In the culmination phase, clients contribute their
learned features to a central aggregation point, culminating in constructing a
Federated Global Model. By virtue of collaborative learning, this model is
anticipated to surpass the predictive capabilities of any single client’s
model, thereby enhancing the efficacy of real-time image classification.

Data Decentralization and Classification: In the initial phase, the framework
undertakes the task of data decentralization, where MRI brain tumor data,
encompassing four distinct classes—Gliomas, Meningiomas, Pituitary tumors, and
non-tumorous images—are distributed among four clients (num_clients=4). This
stratification ensures that each client holds a representative subset of the
entire data spectrum, facilitating the diversity and richness required for a
comprehensive learning process.

Ensemble Model-Based MRI Analysis on a Local Server: The second phase is
dedicated to the local training of models on each client’s server. Where an
ensemble classification approach utilizing $M$ distinct classifiers, denoted
as $h_{m}(X)$, for dataset $X$. The final classification, $C_{X}$, is
determined by a majority vote, or mode, of the outcomes from all classifiers:

$C_{X}=\text{mode}{h_{1}(X),\dots,h_{M}(X)}.$ (1)

Weights $w_{j}$ are assigned to each classifier to influence their
contribution to the final decision, leading to a weighted majority vote
mechanism:

$C_{X}=\arg\max_{i}\sum_{j=1}^{M}w_{j}\mathbb{I}(h_{j}(X)=i),$ (2)

where $\mathbb{I}$ is the indicator function. This method effectively combines
multiple predictions, leveraging the strengths of diverse classifiers. The
ensemble approach’s ability to mitigate the risk of overfitting, thereby
ensuring robust feature extraction and model generalization.

Federated Aggregation and Global Model Development: In the final phase, the
individually trained models from each client undergo a process of federated
aggregation. Here, the unique feature extraction insights gleaned from each
client’s ensemble model are shared with a central aggregation point. This
collaborative effort forms a Federated Global Model, which amalgamates
distributed intelligence into a singular, potent model. Subsequently, this
enhanced global model is redistributed to each client-server, enabling the
prediction of real-time images with augmented accuracy and reliability.

Algorithm 2 introduces a groundbreaking approach to medical image analysis by
combining multiple DL models with FL, enhancing diagnostic accuracy while
preserving privacy and integrity.

Data: MRI_Brain_tumor_data, num_clients=4

Result: Enhanced prediction accuracy for real-time images

1

2

3Function  _Main(__)_ :

4 for  _client_id to num_clients_ do

5 local_data $\leftarrow$ DistributeData(MRI_Brain_tumor_data, client_id)

6 local_model $\leftarrow$ TrainClientModel(_local_data, client_id_)

7 SendToServer(local_model)

8

9 global_model $\leftarrow$ AggregateModels(__)

10 for  _client_id to num_clients_ do

11 MakePrediction(_global_model, client_id_)

12

13

14 Function  _TrainClientModel(_data, client_id_)_:

15 architectures $\leftarrow$ [MobileNetV2, VGG16, VGG19]

16 ensemble_predictions $\leftarrow$ [ ]

17 for  _architecture in architectures_ do

18 model $\leftarrow$ InitializeModel(architecture)

19 trained_model $\leftarrow$ TrainModel(model, data)

20 ensemble_predictions.append(trained_model)

21

22 return CombineModels(ensemble_predictions, method=’voting’)

23

24Function  _AggregateModels(__)_ :

25 aggregated_models $\leftarrow$ ReceiveFromClients()

26 weighted_models $\leftarrow$ CalculateWeights(aggregated_models)

27 global_model $\leftarrow$ WeightedSum(weighted_models)

28 return global_model

29

30

31Function  _MakePrediction(_model, client_id_)_:

32 test_data $\leftarrow$ LoadClientTestData(client_id)

33 prediction $\leftarrow$ model.Predict(test_data)

34 return prediction

35

36

Algorithm 1 EBFL Framework for MRI Analysis

###  III-B Mathematical Formulation for the Proposed Maximum Voting Classifier
Techniques incorporating federated learning

Given a set of $N$ classes for brain tumor classification, let
$\mathcal{C}=\\{c_{1},c_{2},\ldots,c_{N}\\}$ represent the class labels for
glioma, meningioma, pituitary, and no-tumor. For a dataset $X$ with $m$ MRI
images, $X=\\{x_{1},x_{2},\ldots,x_{m}\\}$, and a set of $K$ base deep
learning models $\mathcal{M}=\\{M_{1},M_{2},\ldots,M_{K}\\}$, the ensemble-
enhanced federated learning framework operates as follows:

Data Distribution: The dataset $X$ is partitioned into $P$ subsets
corresponding to the number of clients, such that $X_{p}\subset X$ for
$p=1,\ldots,P$, where $P$ is the total number of clients.

Local Model Training: Each client $p$ trains a local ensemble model $E_{p}$
using its subset of data $X_{p}$, where $E_{p}$ is a combination of models
$\mathcal{M}$ trained on $X_{p}$. The ensemble model at client $p$ is defined
as [40]:

$E_{p}(x)=\text{mode}\\{M_{1}(x),M_{2}(x),\ldots,M_{K}(x)\\}\quad\text{for}\quad
x\in X_{p}$ (3)

Global Model Aggregation: After training, clients contribute their model’s
knowledge to a central server. The server aggregates these updates to form a
global model $G$, which is then distributed back to the clients for prediction
on new data.

Prediction with the Global Model: The global model $G$ is used to predict the
class label $y$ for a new MRI image $x$ as follows:

$y=G(x)=\text{mode}\\{E_{1}(x),E_{2}(x),\ldots,E_{P}(x)\\}$ (4)

The ensemble model’s prediction $y$ for an input instance $x$ is the class
label that receives the majority vote from the predictions of the local
ensemble models $E_{p}$.

For the loss function, the multiclass cross-entropy loss for a dataset with
$m$ samples and $N$ classes is used to optimize the ensemble model and is
defined as [40]:

$L=-\frac{1}{m}\sum_{i=1}^{m}\sum_{j=1}^{N}y_{ij}\log(\hat{y}_{ij})$ (5)

where $y_{ij}$ is the binary indicator (0 or 1) if class label $j$ is the
correct classification for observation $i$, and $\hat{y}_{ij}$ is the
predicted probability of observation $i$ being of class $j$.

Algorithm 2 captures the essence of using multiple classifiers in a federated
learning setting to address the challenge of brain tumor classification,
emphasizing preserving data privacy and utilizing the strengths of ensemble
learning.

Data: MRI_Brain_tumor_data, num_clients=4

Result: Enhanced prediction accuracy for real-time images

1

2

3Function  _Main(__)_ :

4 for  _client_id to num_clients_ do

5 local_data $\leftarrow$ DistributeData(MRI_Brain_tumor_data, client_id)

6 local_model $\leftarrow$ TrainClientModel(_local_data, client_id_)

7 SendToServer(local_model)

8

9 global_model $\leftarrow$ AggregateModels(__)

10 for  _client_id to num_clients_ do

11 MakePrediction(_global_model, client_id_)

12

13

14

15 Function  _TrainClientModel(_data, client_id_)_:

16 architectures $\leftarrow$ [MobileNetV2, VGG16, VGG19]

17 ensemble_predictions $\leftarrow$ [ ]

18 for  _architecture in architectures_ do

19 model $\leftarrow$ InitializeModel(architecture)

20 trained_model $\leftarrow$ TrainModel(model, data)

21 ensemble_predictions.append(trained_model)

22

23 return CombineModels(ensemble_predictions, method=’voting’)

24

25

26Function  _AggregateModels(__)_ :

27 aggregated_models $\leftarrow$ ReceiveFromClients()

28 weighted_models $\leftarrow$ CalculateWeights(aggregated_models)

29 global_model $\leftarrow$ WeightedSum(weighted_models)

30 return global_model

31

32

33Function  _MakePrediction(_model, client_id_)_:

34 test_data $\leftarrow$ LoadClientTestData(client_id)

35 prediction $\leftarrow$ model.Predict(test_data)

36 return prediction

37

38

Algorithm 2 EBFL Framework for MRI Analysis

###  III-C Multi-Class Brain Tumor Classification Using Federated Learning
mathematical explanation

Each client $C_{i}$ trains a local model $M_{i}$ on their own dataset $D_{i}$.
The objective is to minimize the local loss function $L_{i}$, which is a
function of the model parameters $\theta_{i}$ and the data $D_{i}$:

$L_{i}(\theta_{i})=\frac{1}{|D_{i}|}\sum_{(x,y)\in
D_{i}}l(M_{i}(x;\theta_{i}),y)$ (6)

where $l$ is a loss function like cross-entropy for classification tasks, $x$
is the input data, and $y$ is the label.

### Local Model Evaluation

After training, each client evaluates the performance of their model on a
validation set to calculate metrics such as accuracy, precision, recall, and
F1-score.

### Model Aggregation (Federated Averaging)

The server aggregates the parameters $\theta_{i}$ from all clients to update
the global model $M_{g}$ using federated averaging. The new global parameters
$\theta_{g}$ are computed as the weighted average of local parameters:

$\theta_{g}=\frac{\sum_{i=1}^{N}w_{i}\theta_{i}}{\sum_{i=1}^{N}w_{i}}$ (7)

where $w_{i}$ is the weight for client $i$, often chosen as the number of
samples in $D_{i}$, and $N$ is the total number of clients.

##  IV  Experimental Analysis

###  IV-A Dataset Description & experimental result

Our research employs a comprehensive dataset of 6,771 brain MRI scans sourced
from the Kaggle data repository [41] to classify brain tumors. This dataset
encompasses four distinct categories: glioma (1,621 images), characterized by
aggressive tumors arising from the glial cells; meningioma (1,645 images),
typically slower-growing tumors originating from the meninges; pituitary
tumors (1,775 images), often benign growths in the pituitary gland; and non-
tumorous scans (2,000 images), serving as a control group. The dataset is
split into 80% for training and 20% for testing, allowing for comprehensive
analysis and validation of the classification performance. Table 1 shows the
full classifier training and validation results of five different approaches.

TABLE I: Comprehensive Assessment of Model Performance

Algorithms |  | Precision
---
(%)
| Recall
---
(%)
| F1-
---
Score
(%)
| Training
---
Accuracy
(%)
| Training
---
loss (%)
| Validation
---
Accuracy
(%)
| Validation
---
loss
(%)
Global Model (FL) | 94 | 91 | 92 | 95 | 0.005 | 94 | 0.06
MobileNetV2 | 94.19 | 92 | 94 | 94 | 0.06 | 96 | 0.06
VGG16 | 90.62 | 91 | 90 | 92 | 0.17 | 91 | 0.20
VGG19 | 94.15 | 95 | 94 | 96 | 0.11 | 94 | 0.13
| Ensemble
---
Model
96 | 97 | 96 | 97 | 0.03 | 96 | 0.04

###  IV-B Evaluating DL Models in Brain Tumor Classification: A Comparative
Analysis

Our study leverages confusion matrices to assess the performance of deep
learning (DL) models in brain tumor classification, including base models
MobileNetV2, VGG16, VGG19, FL aggregated global model, and an ensemble model.
The confusion matrix Figure 2 and training/validation loss and accuracy
metrics Figure 3 are pivotal in evaluating these models. MobileNetV2
demonstrates high precision (94.19%) and recall (92%), with an F1-score of
94%, indicating effective tumor classification and a balance of sensitivity
and specificity. VGG16, though slightly less precise at 90.62%, maintains a
91% recall and an F1-score of 90%, suggesting reliable detection capabilities
with some scope for improvement in reducing training and validation losses.
VGG19 shows notable recall (95%) and precision (94.15%), achieving an F1-score
of 94%, highlighting its sensitivity and low rate of false negatives. Notably,
the ensemble model surpasses individual models with a precision of 96%, recall
of 97%, and an F1-score of 96%, demonstrating its superiority in reducing
variance and bias and enhancing prediction accuracy for brain tumor
classifications.

Figure 2: Ensemble & Base models confusion matrices for Brain tumor multi-
class classification ) Figure 3: Ensemble & Base models Training and
Validation loss and accuracy for identifying brain tumor

###  IV-C Global model (FL) evaluation matrices analysis

The FL approach has culminated in a global model Figure 4 that showcases
commendable precision and recall, with scores of 93.27% and 91%, respectively.
The model achieves an F1-score of 92%, reflecting its effective learning from
distributed data sources while maintaining privacy and data security. The low
training and validation losses indicate a well-fitted model with good
generalization capabilities on unseen data.

Figure 4: Federated Learning model confusion matrix for brain tumor
classification at the local server )

This section delves into the evaluation of the global model’s performance,
focusing on the analysis of the F1 score and validation loss Figure 4 trends
to gauge the efficacy and reliability of the model over the validation
dataset.

##  V Conclusion

Our study introduces the innovative EBFL framework, seamlessly integrating
Federated Learning (FL) with a custom ensemble of advanced deep learning
architectures. This approach significantly enhances feature extraction
capabilities, surpassing single-model systems’ performance while reducing
bias. The EBFL framework prioritizes patient privacy by focusing on feature
extraction rather than direct data sharing. Our results demonstrate the
model’s exceptional accuracy, with the Global model achieving a 94% precision
rate and the Ensemble model further elevating this to 96%. However, the
framework’s main challenge lies in its computational complexity, presenting a
valuable area for future research to explore and optimize.

## Acknowledgment

This material is based on work supported by the National Science Foundation
Award Numbers 2205773 and 2219658.

## References

  * [1] M. Arif, F. Ajesh, S. Shamsudheen, O. Geman, D. Izdrui, D. Vicoveanu, et al., “Brain tumor detection and classification by mri using biologically inspired orthogonal wavelet transform and deep learning techniques,” Journal of Healthcare Engineering, vol. 2022, 2022.
  * [2] A. Ari and D. Hanbay, “Deep learning based brain tumor classification and detection system,” Turkish Journal of Electrical Engineering and Computer Sciences, vol. 26, no. 5, pp. 2275–2286, 2018.
  * [3] R. Ranjbarzadeh, A. Caputo, E. B. Tirkolaee, S. J. Ghoushchi, and M. Bendechache, “Brain tumor segmentation of mri images: A comprehensive review on the application of artificial intelligence tools,” Computers in biology and medicine, vol. 152, p. 106405, 2023.
  * [4] A. Amin, K. Hasan, S. Zein-Sabatto, D. Chimba, I. Ahmed, and T. Islam, “An explainable ai framework for artificial intelligence of medical things,” arXiv preprint arXiv:2403.04130, 2024.
  * [5] S. Maqsood, R. Damaševičius, and R. Maskeliūnas, “Multi-modal brain tumor detection using deep neural network and multiclass svm,” Medicina, vol. 58, no. 8, p. 1090, 2022.
  * [6] S. Saeedi, S. Rezayi, H. Keshavarz, and S. R. Niakan Kalhori, “Mri-based brain tumor detection using convolutional deep learning methods and chosen machine learning techniques,” BMC Medical Informatics and Decision Making, vol. 23, no. 1, p. 16, 2023.
  * [7] P. Wu, Z. Wang, B. Zheng, H. Li, F. E. Alsaadi, and N. Zeng, “Aggn: Attention-based glioma grading network with multi-scale feature extraction and multi-modal information fusion,” Computers in biology and medicine, vol. 152, p. 106457, 2023.
  * [8] F. Garcea, A. Serra, F. Lamberti, and L. Morra, “Data augmentation for medical imaging: A systematic literature review,” Computers in Biology and Medicine, vol. 152, p. 106391, 2023.
  * [9] A. Kebaili, J. Lapuyade-Lahorgue, and S. Ruan, “Deep learning approaches for data augmentation in medical imaging: A review,” Journal of Imaging, vol. 9, no. 4, p. 81, 2023.
  * [10] A. Biswas, P. Bhattacharya, S. P. Maity, and R. Banik, “Data augmentation for improved brain tumor segmentation,” IETE Journal of Research, vol. 69, no. 5, pp. 2772–2782, 2023.
  * [11] K. Alomar, H. I. Aysel, and X. Cai, “Data augmentation in classification and segmentation: A survey and new strategies,” Journal of Imaging, vol. 9, no. 2, p. 46, 2023.
  * [12] A. Amin, H. Ma, M. S. Hossain, N. A. Roni, E. Haque, S. Asaduzzaman, R. Abedin, A. B. Ekram, and R. F. Akter, “Industrial product defect detection using custom u-net,” in 2022 25th International Conference on Computer and Information Technology (ICCIT), pp. 442–447, IEEE, 2022.
  * [13] L. Fezai, T. Urruty, P. Bourdon, C. Fernandez-Maloigne, and A. D. N. Initiative, “Deep anonymization of medical imaging,” Multimedia Tools and Applications, vol. 82, no. 6, pp. 9533–9547, 2023.
  * [14] M. U. Akbar, W. Wang, and A. Eklund, “Beware of diffusion models for synthesizing medical images–a comparison with gans in terms of memorizing brain tumor images,” arXiv preprint arXiv:2305.07644, 2023.
  * [15] J. Weigand-Whittier, M. Sedykh, K. Herz, J. Coll-Font, A. N. Foster, E. R. Gerstner, C. Nguyen, M. Zaiss, C. T. Farrar, and O. Perlman, “Accelerated and quantitative three-dimensional molecular mri using a generative adversarial network,” Magnetic Resonance in Medicine, vol. 89, no. 5, pp. 1901–1914, 2023.
  * [16] G. A. Kaissis, M. R. Makowski, D. Rückert, and R. F. Braren, “Secure, privacy-preserving and federated machine learning in medical imaging,” Nature Machine Intelligence, vol. 2, no. 6, pp. 305–311, 2020.
  * [17] M. J. Willemink, W. A. Koszek, C. Hardell, J. Wu, D. Fleischmann, H. Harvey, L. R. Folio, R. M. Summers, D. L. Rubin, and M. P. Lungren, “Preparing medical imaging data for machine learning,” Radiology, vol. 295, no. 1, pp. 4–15, 2020.
  * [18] S. Goyal and R. Singh, “Detection and classification of lung diseases for pneumonia and covid-19 using machine and deep learning techniques,” Journal of Ambient Intelligence and Humanized Computing, vol. 14, no. 4, pp. 3239–3259, 2023.
  * [19] N. Nasir, A. Kansal, F. Barneih, O. Al-Shaltone, T. Bonny, M. Al-Shabi, and A. Al Shammaa, “Multi-modal image classification of covid-19 cases using computed tomography and x-rays scans,” Intelligent Systems with Applications, vol. 17, p. 200160, 2023.
  * [20] Z. Ullah, M. Usman, S. Latif, and J. Gwak, “Densely attention mechanism based network for covid-19 detection in chest x-rays,” Scientific Reports, vol. 13, no. 1, p. 261, 2023.
  * [21] S. Al Jannat, A. Amin, M. S. Hossain, E. Hossain, E. Haque, and N. A. Roni, “Adpt: An automated disease prognosis tool towards classifying medical disease using hybrid architecture of deep learning paradigm,” in 2022 25th International Conference on Computer and Information Technology (ICCIT), pp. 903–908, IEEE, 2022.
  * [22] N. Sollmann, M. T. Löffler, S. Kronthaler, C. Böhm, M. Dieckmeyer, S. Ruschke, J. S. Kirschke, J. Carballido-Gamio, D. C. Karampinos, R. Krug, et al., “Mri-based quantitative osteoporosis imaging at the spine and femur,” Journal of Magnetic Resonance Imaging, vol. 54, no. 1, pp. 12–35, 2021.
  * [23] T. Jentzsch, D. W. Cadotte, J. R. Wilson, F. Jiang, J. H. Badhiwala, M. A. Akbar, B. Rocos, R. G. Grossman, B. Aarabi, J. S. Harrop, et al., “Spinal cord signal change on magnetic resonance imaging may predict worse clinical in-and outpatient outcomes in patients with spinal cord injury: a prospective multicenter study in 459 patients,” Journal of clinical medicine, vol. 10, no. 20, p. 4778, 2021.
  * [24] F. R. Santiago, A. J. L. Ramos-Bossini, Y. X. J. Wáng, and D. L. Zúñiga, “The role of radiography in the study of spinal disorders,” Quantitative Imaging in Medicine and Surgery, vol. 10, no. 12, p. 2322, 2020.
  * [25] C. De Vente, P. Vos, M. Hosseinzadeh, J. Pluim, and M. Veta, “Deep learning regression for prostate cancer detection and grading in bi-parametric mri,” IEEE Transactions on Biomedical Engineering, vol. 68, no. 2, pp. 374–383, 2020.
  * [26] Q. Wu, S. Wang, S. Zhang, M. Wang, Y. Ding, J. Fang, W. Qian, Z. Liu, K. Sun, Y. Jin, et al., “Development of a deep learning model to identify lymph node metastasis on magnetic resonance imaging in patients with cervical cancer,” JAMA network open, vol. 3, no. 7, pp. e2011625–e2011625, 2020.
  * [27] B. Abhisheka, S. K. Biswas, and B. Purkayastha, “A comprehensive review on breast cancer detection, classification and segmentation using deep learning,” Archives of Computational Methods in Engineering, vol. 30, no. 8, pp. 5023–5052, 2023.
  * [28] B. Fritz, H. Y. Paul, R. Kijowski, and J. Fritz, “Radiomics and deep learning for disease detection in musculoskeletal radiology: an overview of novel mri-and ct-based approaches,” Investigative radiology, vol. 58, no. 1, pp. 3–13, 2023.
  * [29] N. Kamiya, “Deep learning technique for musculoskeletal analysis,” Deep Learning in Medical Image Analysis: Challenges and Applications, pp. 165–176, 2020.
  * [30] R. Hemke, C. G. Buckless, A. Tsao, B. Wang, and M. Torriani, “Deep learning for automated segmentation of pelvic muscles, fat, and bone from ct studies for body composition assessment,” Skeletal radiology, vol. 49, pp. 387–395, 2020.
  * [31] L. Zhang, M. Wang, M. Liu, and D. Zhang, “A survey on deep learning for neuroimaging-based brain disorder analysis,” Frontiers in neuroscience, vol. 14, p. 779, 2020.
  * [32] N. Yamanakkanavar, J. Y. Choi, and B. Lee, “Mri segmentation and classification of human brain using deep learning for diagnosis of alzheimer’s disease: a survey,” Sensors, vol. 20, no. 11, p. 3243, 2020.
  * [33] A. Amin, K. Hasan, S. Zein-Sabatto, D. Chimba, I. Ahmed, and T. Islam, “An explainable ai framework for artificial intelligence of medical things,” in IEEE Global Communications conference 2023, pp. 1–6, IEEE, 2023, in press.
  * [34] A. Raza, H. Ayub, J. A. Khan, I. Ahmad, A. S. Salama, Y. I. Daradkeh, D. Javeed, A. Ur Rehman, and H. Hamam, “A hybrid deep learning-based approach for brain tumor classification,” Electronics, vol. 11, no. 7, p. 1146, 2022.
  * [35] W. Ayadi, W. Elhamzi, I. Charfi, and M. Atri, “Deep cnn for brain tumor classification,” Neural processing letters, vol. 53, pp. 671–700, 2021.
  * [36] J. D. Bodapati, N. S. Shaik, V. Naralasetti, and N. B. Mundukur, “Joint training of two-channel deep neural network for brain tumor classification,” Signal, Image and Video Processing, vol. 15, no. 4, pp. 753–760, 2021.
  * [37] J. Kang, Z. Ullah, and J. Gwak, “Mri-based brain tumor classification using ensemble of deep features and machine learning classifiers,” Sensors, vol. 21, no. 6, p. 2222, 2021.
  * [38] A. Shafi, M. B. Rahman, T. Anwar, R. S. Halder, and H. E. Kays, “Classification of brain tumors and auto-immune disease using ensemble learning,” Informatics in Medicine Unlocked, vol. 24, p. 100608, 2021.
  * [39] S. Nazir and M. Kaleem, “Federated learning for medical image analysis with deep neural networks,” Diagnostics, vol. 13, no. 9, p. 1532, 2023.
  * [40] G. M. James, “Majority vote classifiers: theory and applications,” 1998.
  * [41] M. NICKPARVAR, “Brain tumor mri dataset,” 2021.  Accessed: 2023-04-08.