\section{Introduction}
With the success of deep neural networks (NNs) in recent years, there has been an increasing trend to introduce machine learning based approaches into many domains -- including safety-critical areas,
such as airborne collision avoidance~\cite{julian2016policy} and autonomous cars~\cite{bojarski2016end}.
Therefore,
there is a growing interest in verification of NNs.
This spurred research on new verification methods~\cite{katz2017reluplex,Singh2018,wang2018efficient,tran2019star,Buning2020,paulsen2020neurodiff,narodytska2018verifying}.
The literature can broadly be classified into adversarial robustness verification (e.g. the work by Singh et al.~\cite{Singh2018}), functional property verification (e.g. the work by Katz et al.~\cite{katz2017reluplex} on the verification of ACAS Xu) and equivalence verification (e.g. Kleine Büning et al.~\cite{Buning2020} and Paulsen et al.~\cite{paulsen2020neurodiff}).
The main application of equivalence verification
is in the space of NN compression.
As NNs grow ever larger and computing becomes ever more ubiquitous,
resource restrictions require to compress large NNs into smaller models.
Cheng et al.~\cite{cheng2020survey} give an extensive survey of such compression techniques.
Furthermore, equivalence verification can be deployed to examine the influence of certain NN-based pre-processing steps (cf.~\cite{narodytska2018verifying}) or in cases where performing multiple verification tasks on a large NNs would be too expensive (cf.~\cite{Buning2020}).
In what follows, we refer to the current method of Paulsen et al.~\cite{paulsen2020reludiff} by \textsc{ReluDiff} and that of Kleine Büning et al.~\cite{Buning2020} by \textsc{MilpEquiv}.

One approach which has previously been shown to yield good results for adversarial and functional verification on NNs is Geometric Path Enumeration (GPE)~\cite{tran2019star,bak2020improved}.
However, this algorithm was initially devised as an approach operating on a single NN.
In this work we extend GPE to a setting with multiple NNs and implement its extension for the problem of equivalence verification.
We explore which (sometimes previously used) optimizations yield good results when applied to the equivalence problem.
While our work in this paper is specific to the problem of equivalence verification, the extended GPE algorithm can also be used for other verification tasks involving multiple NNs.

\subsection{Contribution}
In this work, we focus on the problem of equivalence verification for (potentially) structurally differing NNs.
Our contributions are as follows:
\begin{itemize}
    \item[(C1)] We prove that the $\epsilon$-equivalence problem for NNs is coNP-com\-plete.

    \item[(C2)]
    We extend the GPE algorithm (Tran et al.~\cite{tran2019star}) to a setting with multiple NNs and apply it to the equivalence verification problem.

    \item[(C3)] We evaluate several optimizations for this setting which increase efficiency on practical problems.

    \item[(C4)] We perform a comparative evaluation of our algorithm (on ACAS Xu and modified MNIST benchmarks) and show that it outperforms \textsc{MilpEquiv} in four of our five $\epsilon$-equivalence benchmarks and in counterexample finding.
\end{itemize}

\subsection{Overview}
The structure of our paper is as follows:
In \Cref{sec:prelim-related}, we present related work in the field of NN equivalence and introduce the basic notions of GPE relevant to our paper.
Before explaining the idea behind our extension of GPE to multiple NNs in \Cref{sec:extending}, we show that the $\epsilon$-equivalence problem is coNP-complete in \Cref{sec:np-complete}.
Starting from a naive algorithm, we then evaluate optimizations to enable efficient equivalence verification using GPE in \Cref{sec:optimizing}.
We further explore these
optimizations, in particular the question of good refinement heuristics, in \Cref{sec:tree-explore}.
Finally, in \Cref{sec:eval}, we evaluate our algorithm and show advantages and disadvantages to the current state of the art represented by \textsc{MilpEquiv}
\cite{Buning2020}.

\section{Preliminaries \& Related Work}
\label{sec:prelim-related}
\input{sections2/related_work.tex}

\section{NN Equivalence and NP completeness}
\label{sec:np-complete}
\input{sections2/np_complete.tex}

\section{Extending GPE to multiple NNs}
\label{sec:extending}
\input{sections2/extending_gpe.tex}

\section{Optimizing GPE for two NNs}
\label{sec:optimizing}
\input{sections2/optimizing_gpe.tex}

\section{The Branch Tree Exploration Problem}
\label{sec:tree-explore}
\input{sections2/branch_tree_exploration2.tex}

\section{Experimental Evaluation}
\input{sections2/experiments.tex}

\section{Conclusion and Future Work}
\label{sec:conclusion}
We proposed an approach extending Geometric Path Enumeration~\cite{tran2019star} to multiple NNs.
Employing this method, we presented an equivalence verification algorithm which was
optimized by four techniques: Zonotope propagation, zonotope over-approximation, LP approximation, and refinement heuristics.
Our evaluation shows that the optimizations increase the approach's efficiency and that it can verify equivalence of NNs which were not verifiable by \textsc{MilpEquiv}~\cite{Buning2020} and \textsc{ReluDiff}~\cite{paulsen2020reludiff}.
Our approach significantly outperforms the state of the art in counterexample finding by solving 890 instances in comparison to 315 instances solved by \textsc{MilpEquiv}.
In addition, we proved the coNP-completeness of the $\epsilon$-equivalence problem, and presented a formal way of reasoning about refinement heuristics in the context of GPE.

In terms of efficiency, one could further explore possible refinement heuristics and consider parallelized (possibly GPU based) implementations.
Moreover, while GPE can increase the confidence in NNs, the role of numerical stability for the verification approach has to be further investigated.
Furthermore, an integration of MILP constraints into GPE propagation could be explored resulting in an algorithm inbetween \textsc{NNEquiv} and \textsc{MilpEquiv}.
Additionally, we see a need for a larger body of equivalence benchmarks which allows the conclusive evaluation of equivalence verification algorithms.

\bibliographystyle{IEEEtran}

\section{Introduction}
\textbf{Storyline}
\begin{itemize}
    \item We are concerned with equivalence verification
    \item Equivalence Verification is useful
    \item Unfortunately disproving equivalence is NP complete
    \item Previous approaches differ in structurally similar/structurally different
    \item We are interested in structurally different
    \item We extend GPE to a multi-network usecase
    \item Here are the optimizations we came up with
    \item Those are the open questions (in particular optimal? heuristics)
    \item While we cannot beat the state of the art globally we can prove certain things faster than \textcite{Buning2020}
\end{itemize}
In recent years neural networks have been increasingly introduced to safety critical domains...
\begin{itemize}
    \item Why do we need equivalence verification?
    \item Distinguish structurally similar to structurally different networks
    \item For the one kind we have the approach presented by \textcite{paulsen2020neurodiff}. For the other kind we have the approach presented by \textcite{Buning2020}.
    \item However there are still scalability issues, which begs the question whether approaches previously successful in single neural network verification such as GPE\cite{tran2019star,bak2020improved} can equally be used for equivalence verification
    \item Single Neural Network verification is NP complete \textcite{katz2017reluplex}. Is that the case for epsilon/top1 equivalence verification?
\end{itemize}
\subsection{Contribution}
\todo{Unfinished, needs reworking}
\begin{itemize}
    \item We extend the basic GPE algorithm to the be applicable for equivalence verification
    \item We introduce the necessary optimizations which allow GPE to work on this problem set. Mention increased exponential growth to larger number of ReLU nodes
    \item We experimentally compare our algorithm to the current State of the Art for equivalence verification for structurally differing neural networks\cite{Buning2020}, analyze currently remaining deficiencies and bottlenecks of our algorithm and present possible pathways for further optimizations.
\end{itemize}
\subsection{Overview}

\section{Preliminaries & Related Work}
\input{sections/related_work.tex}

\section{Extending GPE to multiple networks}
\input{sections/extending_gpe.tex}

\section{Optimizing GPE for two networks}
\input{sections/optimizing_gpe.tex}

\section{The Branch Tree Exploration Problem}
\input{sections/branch_tree_exploration.tex}
\textbf{Notes:}
\begin{itemize}
    \item Right, and what about when to overapproximate?
    \item Yay another problem! Zonotope overapproximation produces more overapproximation, produces more error. We need to change that. Well, how about we use CEGAR? No that has already been shown to be bad by \textcite{10.1007/978-3-030-76384-8_2}, \textbf{but} their approach of having a fixed threshold and mean requires parameters we don't have.
    Let's look at the depth necessary by CEGAR... Oh, look at that! We can see something like a pattern here, let's apply some naive heuristics on that...
    \item Maybe formalize the issue we are currently having
    \item We have a tree of unknown structure
    \item At every node we can decide to try an overapproximation at cost $c_o$.
    \item If successful we are done with this branch, if not, or we decide against we need to descend at cost $c_d$
    \item We know somethings about the tree such as that the overapproximation success function is monotone or that the probability of success increases with depth.
    \item more?
    \item What's a good strategy here?
\end{itemize}

\section{Experimental Evaluation}
\input{sections/experiments.tex}
Benchmarks:

\todo{Clean up Bibliography, both works by Paulsen are published now...}
\printbibliography

\subsection{NEW}
For given networks $\mathcal{R},\mathcal{T}$ and an input space $I$ we define a tree $T$ consisting of all branches the exact Zonotope GPE algorithm explores for its equivalence verification.
We denote $T=\left(V,l,r\right)$ as the tree where $V$ corresponds to the inner nodes and leafs of the tree and $l,r:V\to V\cupplus\left\{\bot\right\}$ denote left and right successors of some node in $V$ such that the structure induced by $V,l$ and $r$ produces a binary tree with some root $S$.
Further, $V$ can be partitioned into $V=N\cupplus R \cupplus W$ where $N$ corresponds to the inner nodes with $\forall v \in N. l\left(v\right)\neq \bot \neq r\left(v\right) \land l\left(v\right)\neq r\left(v\right)$ and $R$ (resp. $W$) corresponds to the feasible (resp. infeasible) leafs $\forall v \in R \cupplus W. l\left(v\right) = r\left(v\right) = \bot$.
An exact execution of GPE implicitly produces a set of paths of the form $p \in N^*\times\left(R \cupplus W\right)$ that are explored sequentially.
This can be written as $P=\left\{p_1, \dotsi, p_{\left|R \cupplus W\right|}\right\}$.
Note that the number of paths produced is fixed to the number of feasible or infeasible leafs.
Additionally it is worth noting that every leaf can be mapped to a subset of the input space which is being verified on this leaf, i.e. we can define a function $A: R \cupplus W \to 2^I$ which maps a leaf to a subset of $I$ (in particular $\forall w \in W. A\left(w\right)=\emptyset$).
We can then lift this function to all nodes in $V$ by defining $A$ for all nodes in $N$ through $A\left(v\right) = A\left(l\left(v\right)\right) \cup A\left(r\left(v\right)\right)$.
In this case the root node of the tree $S \in N$ must fulfill $I \subseteq A\left(S\right)$.

The time taken to verify using exact Zonotope GPE is bounded through:
\[
\mathcal{O}\left( \left|P\right|*t_{\text{check}} + \sum_{p \in P} \left|p\right|*t_{\text{descend}} \right)
\]
Essentially, there is only one way in which we can speed up computation time for GPE:
By finding a tree $T'$ which has a smaller number of leafs $\left|R'\cupplus W'\right|$ or less deep paths $\left|p\right|$, but is also equivalent to $T$ in the sense that the root node $S'$ of $T'$ also has the property $I \subseteq A\left(S\right)$.
However, the problem is further complicated by the property which we want to prove.
Otherwise, we could otherwise simply choose a one leaf tree which covers all of $I$.
For a given equivalence property $\mathcal{E}$ we can thus define a function $\min_{\mathcal{E}}:N^*\times\left(R \cupplus W\right) \to \mathbb{N}$ which returns the number of necessary steps in the given path $p$ for the property $\mathcal{E}$ to be verifiable on $A\left(p_{\min_{\mathcal{E}}\left(p\right)}\right)$.
The exploration of the tree induced by the set of paths
\[
P^*=\left\{ p' \in N^*\times\left(R \cupplus W\right) | \exists p \in P: p'=p_{0 : \min_{\mathcal{E}}\left(p\right)}\right\}
\]
would then be sufficient to prove equivalence ($p_{0:\min_{\mathcal{E}}\left(p\right)}$ denotes the prefix path of $p$ from step $0$ to step $\min_{\mathcal{E}}\left(p\right)$).
While $\left|P\right|$ is the upper bound for the number of paths we have to consider, $\left|P^*\right|$ is the lower bound (in this work we omit the option of reordering the inner nodes of the tree which could, potentially further compress the tree).
Since GPE still needs to process the nodes on which no split is performed (i.e. it needs to compute the overapproximation), we denote the cost of this processing as $t_{\text{overapprox}}$ which must be added to the total cost for any node on the path which is not split. For our model we assume $t_{\text{overapprox}}<t_{\text{descend}}$ since the split needs to replicate the entire Zonotope and its connected state.
We can now take a more formal look at why CEGAR performs particularly bad.
CEGAR begins by computing an overapproximation based on the tree's root node and then tries to prove the property $\mathcal{E}$ at hand. If it fails, it branches on the root node and retries.
This strategy is performed iteratively until each path is successful.
Thus, CEGAR reaches the minimum number $P^*$ of paths to consider.
However, it comes with serious drawbacks: The check function is called significantly more often.
Additionally, for the $i$th refinement $p-i$ of the path's nodes need to be over-approximated.
Through the Gaussian sum we obtain a quadratic number of over-approximations.
This leads to the following lower bound for the algorithm's cost:
\[
\Omega\left(\sum_{p\in P^*} \left|p\right|*(t_{\text{check}}+t_{\text{descend}}) + \left|p\right|^2*t_{\text{overapprox}}\right)
\]
Note, that this is a lower-bound because we omit the cost $t_{\text{overapprox}}$ for the nodes which are over-approximated after we reach the necessary level of over-approximation in $P^*$.
Even when assuming the approximation to be free (i.e. $t_{\text{overapprox}}=0$) this approach becomes less effective than exact propagation once $\sum_{p\in P^*} \left|p\right| > \left|P\right|$.
While there are cases where this approach is effective, this is not always guaranteed to be the case -- especially for larger networks.

The optimal approach would use an oracle which tells the algorithm the number of refinement steps necessary for verifying $\mathcal{E}$ on the given path (i.e. $\min_{\mathcal{E}}$.
Assuming $t_{\text{overapprox}}=0$, the algorithm's cost would reduce to:
\[
\mathcal{O}\left(\left|P^*\right|*t_{\text{check}} + \sum_{p\in P^*} \left|p\right| * t_{\text{descend}}\right)
\]
It follows from the definition of $P^*$ that $\left|P^*\right| \leq \left|P\right|$ and that the average path length of $P^*$ is certain to be shorter or equal to the average length in $P$.
This implies that the oracle based algorithm would be certain to perform better than exact Zonotope GPE.
\todo{We could support this argument by experimentally running the oracle version (step 1: run CEGAR, step 2: run oracle version)}

It is thus important to find a good trade-off between the number of paths to consider and the number of operations necessary per path.
While the exact Zonotope approach yields a large number of paths at a comparatively low cost per path, the CEGAR approach yields the minimal number of paths at a very high cost per path.
We therefore explored a number of heuristics which are trying to estimate the minimum path length.
Again assuming $t_{\text{overapprox}}=0$ it is clear that a too low estimate produces additional checks while a too high estimate produces additional paths which need to be verified -- both either influencing the cost per path or the number of paths considered.
To this end we explored the following options:
\begin{enumerate}
    \item Maximum?
    \item Average*x?
    \item Running Average?
    \item Running Maximum?
    \item Running 95 percentile?
    \item Last Value
    \item Last Value - 1
\end{enumerate}

These strategies are much more difficult to analyze theoretically because they are particularly dependent on the distribution of the encountered paths and would thus only allow a probabilistic analysis for which we would then have to assume a probability distribution on the paths.

Maybe instead just give a rationale why we think certain strategies work better or worse?

\subsection{OLD}
We consider the algorithm in \Cref{algorithm:tree-exploration} which explores a binary tree.
At every node the algorithm has the possibility to either further descend into the tree in a depth first manner, or the algorithm can attempt to use an over-approximation.
This decision is made by the function \texttt{do\_approx}.
If the over-approximation is successful (executed and returned by \texttt{approx\_succ}, the algorithm must not further explore the current branch.
On the other hand, if the over-approximation fails, we must further \texttt{descend} into the children of the current node.
Both operations, descend and over-approximation check, come at a cost and it is thus necessary to find an optimal strategy which reduces the total cost of the tree exploration.
The tree nodes correspond to the neurons which we may either split on or over-approximate.
\todo{Only include this algorithm if it actually turns out to be interesting of course...}
\begin{algorithm}
    \caption{Tree Exploration Algorithm}
    \label{algorithm:tree-exploration}
    \begin{algorithmic}
    \REQUIRE Tree root $R$
    \ENSURE All branches explored or over-approximated
    \STATE $S \leftarrow \langle R \rangle$
    \WHILE{$\neg\texttt{S}.\texttt{empty}$}
        \STATE $n \leftarrow \texttt{S.pop}\left(\right)$
        \IF{$\texttt{do\_approx}\left(n\right)$\&\&$\texttt{approx\_succ}\left(n\right)$}
            \STATE \textbf{continue}
        \ENDIF
        \STATE $l,r \leftarrow \texttt{descend}\left(n\right)$
        \STATE $\texttt{S.pushAll}\left(\langle l,r \rangle\right)$
    \ENDWHILE
    \end{algorithmic}
\end{algorithm}
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/acas-succ-depth.png}
    \caption{ACASXU\_run2a\_1\_1: Descend depth at point of successful equivalence proof for $\epsilon=0.05$ (Running percentile window: 479)}
    \label{fig:success-depth}
\end{figure}
In order to find the optimal positions for an over-approximation, we can execute the previously outlined naive refinement loop which iteratively refines more and more nodes from the start.
We can then track the number of \texttt{descend} operations which were necessary before \texttt{approx\_succ} was successful.
Consider the example of two ACAS Xu networks, a benchmark previously used for equivalence verification by \textcite{paulsen2020neurodiff}.
In \Cref{fig:success-depth} we show the number of necessary splits and some running percentiles for the Zonotopes considered during the equivalence verification of an ACAS Xu network, previously considered for equivalence verification by \textcite{paulsen2020neurodiff}.\todo{Since we talk about how \textcite{10.1007/978-3-030-76384-8_2} only used ACAS, we may want to use an MNIST -- at least additionally}.
The x axis represents the Zonotopes sequentially encountered over the runtime of the algorithm while the y axis describes the number of necessary refinement steps.
As can be seen from the data there is a high variance in the number of splits even locally. However, at the same time we see that the moving percentiles vary over time, thus indicating that certain \textit{phases} of the path enumeration require more splits than others.
This raises the question how the necessary number of refinement steps can be analyzed or, ideally, predicted based on previously encountered data.

Approaches considered so far:
\begin{itemize}
    \item Compute overall maximum. Problem: Does not fall once maximum was reached
    \item Take Last
    \item Take Last - 1
    \item Take Last - 2 (this one seemed worse)
    \item Take 95 percentile of last 300 or so(?)
    \item Taking running max of last 300 or so
\end{itemize}
\subsection{Comparison of Optimizations \& Heuristics}
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/preliminary-comparison.png}
    \caption{Time taken for the 8 equivalent benchmarks considered:...}
    \label{fig:my_label}
\end{figure}
Benchmarks: A1, A2, D, E, F
\begin{itemize}
    \item Naive exact star
    \item Zono
    \item Zono+LP feasibility (i.e. exact)
    \item Zono Overapprox+LP approximation
    \begin{itemize}
        \item CEGAR
        \item Optimal
        \item Last-1
        \item Running Max
    \end{itemize}
\end{itemize}

\subsection{$\epsilon$ tightness}
\todo{Insert plot}
\textbf{NOTE}
Why did NNEquivalence return that bound while I found a counterexample?

\subsection{Counterexample finding}
\textbf{NOTE}
This one is difficult...
I mean I could easily find 10 examples where NNEquiv "instantly" returns a counter-example.
What seems much harder is to find cases where it works a very long time and eventually does find a counter-example.
However logic dictates these cases should exist...
Should we include this anyway?
Should I try to "manipulate" a network to show this behavior?

\subsection{Gurobi vs GLPK?}
I have done \textbf{some} parameter tuning (using gurobi's tune function), but there's certainly more that could be done here...
Use optimal approach from previous section and compare
Benchmarks: A1, A2, D, E, F

\subsection{Comparison to \textcite{Buning2020}}
Top 1:

Depending on input space size we're better or worse...

\subsection{Equivalence on Set Data Structures}
While step functions can, in principle, be reused from previous works on GPE, it becomes necessary to define an equivalence check \texttt{is\_equiv} comparing two set data structures.
To this end we will, first, present how $\epsilon$-equivalence with Chebyshev norm $\left| \cdot \right|_\infty$ can be proven for Zonotopes.
Afterwards we will show how Star Sets can be used to prove a tighter variant of $\epsilon$-equivalence and top-1 equivalence.

\paragraph{$\epsilon$ Equivalence}
In order to prove $\epsilon$ equivalence with the Chebyshev norm, we need to bound the maximum deviation between the two network outputs by $\epsilon$.
That is, given the two output Zonotopes $\Psi^o_1 = \langle c^o_1, G^o_1, P^o_2 \rangle$ and $\Psi^o_2 = \langle c^o_2, G^o_2, P^o_2 \rangle$ we want to find the maximal deviation:
\begin{align*}
    &\max_{\alpha \in P^o_2} \left|\left(c^o_1+G^o_2 \alpha\right) - \left(c^o_2 + G^o_2 \alpha\right)\right|_\infty\\
    = &\max_{\alpha \in P^o_2} \max_{i} \left|\left(c^o_1-c^o_2\right)_i + \left(G^o_1-G^o_2\right)_i \alpha\right|.
\end{align*}
As can be seen by the reformulation above, we can find the maximal deviation over the output by solving optimization problems for each dimension of the \textit{differential Zonotope}
\[
\partial\Psi = \langle \left(c^o_1-c^o_2\right), \left(G^o_1-G^o_2\right), P^o_2 \rangle).
\]
Recalling, that Zonotopes can be optimized with a closed formula, this enables a relatively quick check for the adherence of the desired $\epsilon$ equivalence property.
However, since Zonotopes only approximate the output set, one may need to fall back to the use of Star Sets if equivalence cannot be established using Zonotopes.
In this case, we can reuse the same formula form above to obtain a \textit{differential star set} $\partial\Theta$ which is then optimized using LP solving.

\paragraph{Top 1 Equivalence}
For the case of Top 1 equivalence there are two possible approaches which both rely on propapagated star sets.
On the one hand, we can reuse the encoding by \textcite{Buning2020} using a MILP solver.\todo{Will we try this?}
If, on the other hand, we prefer to use a simplex solver, we can split up the output star sets further thus enabling us to establish Top 1 equivalence.
In this case we solve one optimization problem per output dimesion $j$.
To this end, for an $O$ dimensional output, we generate $O-1$ constraints which ensure that that dimension $j$ is the maximum in $\Theta^o_1$ and add those constraints to our second output $\Theta^o_2$.
We then compute the maximal difference between output dimension $j$ and the other dimensions in the modified version of $\Theta^o_2$.
In this case we must ensure that the maximal difference is always below $0$ to guarantee top 1 equivalence.
Top k equivalence follows trivially by allowing a certain number of dimensions above $0$.

\subsection{Challenges of the naive approach}
While the techniques outlined above, permit a straight forward extension of GPE to multiple networks and thus allow achieving equivalence verification, the approach comes with a number of pitfalls which one must avoid for the approach to succeed.
The most obvious of which is probably the, possibly, increased exponential growth in star sets.
As previously noted, the exact GPE approach based on star sets, splits the star sets on ReLU nodes and thus can yield a high rate of exponential growth.
\textcite{tran2019star} rightly noted that the observed growth usually drastically falls behind the worst case, however the increase in ReLU nodes certainly leads to an increase in necessary splits.
This is particularly the case for ReLU nodes which cut off very similar hyperplanes (such as the two ReLU nodes in a network at the same position with truncated weights).
This can not only double the work, but it may sometimes also lead to precision problems with typically used LP solvers which tend to show problematic behavior when encountering a problem which has two "nearly inversed" hyperplanes as bounds\footnote{In one case the solver would return drastically differing maximum values for the same optimization problem depending on the previous requests or would suddenly deem the problem infeasible.}.
To avoid such numerical problems we thus use 64bit floats by default and always ensure that feasibility is checked at least once by an exact (i.e. rational) LP solver before a branch is declared infeasible.
While these initial improvements help in making GPE for equivalence \textit{possible}, this approach is not yet \textit{scalable}.
Hence, we will devote the next section to various optimizations which we applied to improve the scalability of the approach.
\subsection{Zonotope Propagation}
The first, most obvious, optimization is the use of overapproximating Zonotopes instead of exact star sets.
To this end, we reuse the propagation algorithm previously proposed by \textcite{bak2020improved} which splits the Zonotope on $\relu$ nodes overlapping $0$ and contracts the Zonotope bounds based on the hyper plane that is currently being split.
While this approach reduces the time spent on optimization, it comes at the cost of an increased number of splits. This is, because the Zonotope is not an exact approach. This is particularly worrysome if the Zonotope is propagated through two networks as it may reproduce very similar splits in both networks thus leading to more Zonotopes which need to be considered.
\todo{Maybe we could show some example here. I believe ACAS Xu 1010 had like 75\% splits which were actually infeasible}

\paragraph{Reintroducing LPs for feasibility checks}
However, we can avoid this issue without reintroducing the previous time spent on optimization by carrying a star set with the Zonotope which collects the so far encountered constraints.
This is also useful, because we need to solve LP problems during the top 1 equivalence check which can be done with greater precision if the split constraints are taken into account.
To mitigate this problem we propagate, both, the overapproximate Zonotope and the exact star set.
This is similar to the strategy used by \textcite{bak2020improved}.
However, instead of using the Zonotope as a prefilter, which decides on which ReLU nodes it is even necessary to execute LP optimization, we do not use the star set for optimization during the propagation.
Instead, we only regularly check the feasibility of the constraint set under consideration and abort the further propagation of the datastructure in case the constraints have become unsatisfiable.
In terms of algorithmics, this only executes the first phase of the Simplex algorithm and thus avoids the, sometimes costly, optimization once an initial basis has been found.

\subsection{Zonotope Overapproximation}
Even with these initial optimizations the algorithm is still left with a large number of Zonotopes which have to be evaluated.
We thus have two options to further optimize the algorithm: We can either further reduce the time spent per Zonotope or we can try to reduce the number of Zonotopes which have to be considered.
In order to achieve the second objective, we can over-approximate certain $\relu$ splits through a methodology first presented by \textcite{Singh2018} and later reused by \textcite{10.1007/978-3-030-76384-8_2}: The idea is to introduce an additional dimension to the Zonotope and use it to form a parallelogram around the area covered by the $\relu$ node.
While this decreases the amount of Zonotopes to consider, it also increases the over-approximation error.
In particular, the over-approximation can lead to more $\relu$ nodes supposedly needing to be split.
In \Cref{fig:gpe-splitting} we compare the share of nodes per layer that are ever considered for a split or over-approximation for splitting Zonotope propagation in comparison with overapproximating Zonotope propagation for an ACAS Xu network:
As can clearly be seen in the plot, the overapproximation approach produces a lot more alleged splits than the exact approach.
These \textit{ghost splits} cannot be easily distinguished from actual, necessary splits.
Note that this comparison aggregates the number of splits over all splitted Zonotopes for the exact case -- for a single propagated Zonotope the number of splits is usually even lower.
This can in turn result in even more over-approximation error which \textcite{10.1007/978-3-030-76384-8_2} referred to as \textit{error snowball}.
Since the over-approximation of all $\relu$ nodes can result in a situation where the desired properties are no longer provable due to approximation error, it becomes necessary to find a suitable refinement heuristic which decides what nodes are refined if an equivalence property cannot be established.
We will discuss the search for a suitable refinement strategy in the next section, however there is one more efficiency improvement related to the dimensionality of the LP problems solved.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/split-neurons-acas.png}
    \caption{Share of $\relu$ nodes per network layer which are ever split during GPE. (Overapprox Net1 is hidden behind Overapprox Net2)}
    \label{fig:gpe-splitting}
\end{figure}

\subsection{LP approximation}
With the introduction of over-approximation we encounter an additional problem: As splitting hyperplanes are now no longer dependent on the input variables \textit{only}, but also depend on the dimensions introduced through the over-approximation, it raises the question how to handle this situation in the propagated star set: Should the we equally introduce the star sets dimensionality and thus the dimensionality of the LP problem?
This strategy poses the problem that can increase the computational load on the LP solver.
Instead, we opted to overapproximate the LP problem.
Classically, for an $m$ dimensional Zonotope with initial input dimensionality $I$ we observe a hyperplane cut of the following form:
\begin{align}
\sum_{i=1}^I g_i \alpha_i + \sum_{i=I+1}^m g_i \alpha_i \leq c
\label{eq:exact-lp}
\end{align}
We can now over-approximate this pro inequality by computing $\mu = \min_{\alpha} \sum_{i=I+1}^m g_i$
and constraining the LP problem with the following inequality:
\begin{align}
    \sum_{i=1}^I g_i \alpha_i \leq c - \mu
    \label{eq:approx-lp}
\end{align}
Since any solution for \Cref{eq:exact-lp} implies a solution for \Cref{eq:approx-lp} the second inequality is an over-approximation and can be used to reduce the number of dimensions the LP solver has to handle despite the use of over-approximation on the GPE side.

\subsection{Verifiying Neural Network Equivalence}
There has recently been a line of work which proposed various compression techniques for neural networks (a full review of such the techniques was recently presented by \textcite{cheng2020survey}).
While such techniques have been shown to be useful in practice, most lack a formal proof of their correctness and rely on empirical evidence.
The usage of such techniques thus raises the question how to prove that two neural networks $\mathcal{R}$ and $\mathcal{T}$ or their corresponding mathematical functions $g_{\mathcal{R}}:\mathbb{R}^I\to \mathbb{R}^O, g_{\mathcal{T}}:\mathbb{R}^I\to \mathbb{R}^O$ are equivalent, i.e. that they produce the same results.
To this end, \textcite{Buning2020} examined three possible definitions of equivalence in a predefined subset $X\subseteq \mathbb{R}^I$, two of which we will review here:\todo{Add top k?}
\begin{definition}[$\epsilon$-Equivalence\cite{Buning2020}]
\label{def:epsilon-equiv}
Two networks $\mathcal{R}$ and $\mathcal{T}$ are $\epsilon$-equivalent with respect to a norm $\lVert \cdot \rVert$, if $\left\lVert g_{\mathcal{R}}\left(x\right)-g_{\mathcal{T}}\left(x\right)\right\rVert \leq \epsilon$ for all $x \in X$
\end{definition}
\begin{definition}[One-Hot Equivalence\cite{Buning2020}]
\label{def:one-hot-equiv}
Two networks $\mathcal{R}$ and $\mathcal{T}$ are one-hot equivalent, if $\argmax_i \left(g_{\mathcal{R}}\left(x\right)\right)_i = \argmax_i \left(g_{\mathcal{R}}\left(x\right)\right)_i$ for all $x \in X$
\end{definition}
While Definition \ref{def:epsilon-equiv} is particularly suitable for regression tasks as it can show a very strong form of equivalence, the latter definition represents a more relaxed interpretation of equivalence which can be especially useful for classification tasks.
\textcite{paulsen2020reludiff} refer to the property in Definition \ref{def:epsilon-equiv} as \textit{differential equivalence}.

\textcite{paulsen2020neurodiff} and \textcite{Buning2020} proposed two fundamentally different approaches to neural network equivalence verification with differing application areas.
While \textcite{paulsen2020neurodiff} propose a technique which uses (symbolic) interval propagation on structurally similar neural networks (i.e. networks with the same architecture but differing weights) to prove $\epsilon$ equivalence, \textcite{Buning2020} propose a technique based on Mixed Integer Linear Program solving which encodes (potentially) structurally differing neural networks and the desired equivalence property into an optimization problem.
While the approach by \textcite{paulsen2020neurodiff} was shown to be efficient for cases where a networks weights had been truncated, \textcite{Buning2020} showed that their approach can works well for a number of small networks and all previously mentioned equivalence properties (e.g. two 64 pixel MNIST networks with a total of 84 $\relu$ nodes).
In particular, \textcite{Buning2020} were able to provide a maximal radius around a datapoint for which the property still held.

\subsection{Geometric Path Enumeration}
Geometric Path Enumeration is a methodology originally proposed by \textcite{tran2019star} for verifying safety properties in neural networks.
Given a neural network $\mathcal{N}$, a set of input instances $\mathcal{I}\subseteq\mathbb{R}^I$ and a safety specification $\mathcal{S}\subseteq\mathbb{R}^O$, safety verification concerns with the question whether there exist any instances $i\in\mathcal{I}$ such that $g_{\mathcal{N}}\left(i\right) \notin \mathcal{S}$.
In their work, a safety specification $\mathcal{S}$ is defined as a set of linear constraints on the neural network outputs $\left\{y^{(L)} | C y^{(L)} \leq d\right\}$

\begin{algorithm}
    \caption{High-level path enumeration algorithm for neural networks as described in \cite{bak2020improved}}
    \begin{algorithmic}
    \label{algorithm:gpe-basics:high-level}
    \REQUIRE Input Set $\mathcal{I}$, Unsafe Set $\mathcal{U}$
    \ENSURE Verification result (\texttt{safe} or \texttt{unsafe})
    \STATE $s \leftarrow \left\langle \texttt{layer}:0, \texttt{neuron}: \texttt{None}, \Theta: \texttt{convert}\left(\mathcal{I}\right) \right\rangle$
    \STATE $\texttt{W} \leftarrow \texttt{List}()$ \COMMENT{List of set datastructures to process}
    \STATE $\texttt{W}.\texttt{put}\left(s\right)$

    \WHILE{$\neg\mathcal{W}$.\texttt{empty}()}
        \STATE $s \leftarrow \texttt{W}.\texttt{pop}()$
        \STATE $\texttt{step}\left(s, \texttt{W}\right)$
        \COMMENT{This call may modify \texttt{W} and $s$}
        \IF{$\neg \texttt{finished}\left(s\right)$}
            \STATE $\texttt{W.push}\left(s\right)$
        \ELSIF{$\neg \texttt{verify}\left(s, \mathcal{U}\right)$}
            \RETURN \texttt{unsafe}
        \ENDIF
    \ENDWHILE
    \RETURN \texttt{safe}
    \end{algorithmic}
\end{algorithm}

Instead of pushing single instances through the network and checking whether they satisfy the required safety property, geometric path enumeration feeds an entire set into the neural network and then evaluates whether any parts of the output sets lie outside the safety specification.
The general setup of this idea can be observed in Algorithm \ref{algorithm:gpe-basics:high-level}: The algorithm's input consists of an input set to evaluate and an unsafe set of states that should not be reached from the input set under consideration.
We then convert the input set into an internal data structure storing both a representation of the set and how far the set has been pushed through the network so far.
We add this data structure to a list \texttt{W} of sets to consider.
We then continuously push the data structures inside \texttt{W} through the neural network at hand using the \texttt{step} function until there are either no further sets to consider or an unsafe state has been observed in the network's final layer.

As data structure to be propagated throughout the network, the authors propose to use generalized star sets:
\begin{definition}[Generalized Star Set\cite{tran2019star}]
A generalized star set $\Theta$ is a tuple $\left\langle c, G, P\right\rangle$ where $c\in\mathbb{R}^n$ is the center, $G=\left\{g_1,\cdots,g_m\right\}$ is a set of $m$ vectors in $\mathbb{R}^n$ called basis vectors, and $P:\mathbb{R}^m \to \left\{\top,\bot\right\}$ is a predicate which is (in the works under consideration) restricted to be a conjunction of linear constraints.
The set represented by $\Theta$ is then defined as:
\[
\left\llbracket \Theta \right\rrbracket = \left\{ x \middle\vert x = c + \sum_{i=0}^{m-1} (\alpha_i v_i)\text{ such that } P\left(\alpha_1,\cdots,\alpha_m\right) \right\}
\]
\todo{indices in sum and predicate don't match}
We will sometimes also refer to $G$ as generator and as matrix $G=\left(g_1g_2\cdots g_m\right)$.
\label{def:gpe:star-set}
\end{definition}
Alternatively, it is possible to use Zonotopes which further restrict the type of predicate allowed:
\begin{definition}[Zonotopes\cite{bak2020improved}]
A Zonotope $\Psi$ is a generalized star set $\langle c, G, P \rangle$ with the further restriction that $P$ may only be an interval constraint (i.e. $P$ only enforces a lower and upper bound for each dimension).
\end{definition}

\Cref{algorithm:gpe-basics:high-level} can be completed by providing suitable \texttt{step}, \texttt{finished} and \texttt{verify} functions for the desired data structure as previously done by \textcite{tran2019star}.
Note that minimization/maximization for Generalized Star Sets is much more expensive than optimization for Zonotopes, which can be computed analytically.
While affine transformations can trivially be applied to both data structures, it is less clear how $\relu$ nodes should be handled by the datastructure.
To this end, one can either split the Zonotope into two cases or one can use various over-approximation techniques\cite{Singh2018,10.1007/978-3-030-76384-8_2}.
We call the approach which splits generalized star sets \textit{exact GPE} and other approaches \textit{approximate GPE} since Zonotopes are usually unable to provide an exact output.

\subsection{Experimental Setup}
The benchmark landscape for the task of equivalence verification is still very limited.
Paulsen et al.~\cite{paulsen2020reludiff} proposed a number of benchmark NNs consisting of pairs of NNs differing only in the bit-width of the weights (32 bit vs. 16 bit).
As discussed before, we see this as a restricted use case and are more interested in generic NNs with varying structures and weights.
This is
why we omit a comparison on these NNs where the approach by Paulsen et al.~\cite{paulsen2020neurodiff}  is clearly faster and more precise.
Structurally differing NNs have been previously proposed by Kleine Büning et al.~\cite{Buning2020} who examined 3 NNs of differing layer depths for digit classification on an MNIST~\cite{lecun1998gradient} data set with reduced resolution~\cite{Dua:2019} (8x8 pixels).

In order to
evaluate and compare the approaches we thus proceeded as follows:
First, we decided to look at two types of NNs: Image classification on the 8x8 pixel MNIST data set and NNs used in control systems in the context of an Airborne Collision Avoidance System (ACAS Xu~\cite{julian2016policy}).
Then, based on the original ACAS Xu NNs 1\_1 and 1\_2, we contructed a total of 4 mirror NNs through retraining (ACAS\_1\_1-retrain, ACAS\_1\_2-retrain) and student teacher training~\cite{hinton2015distilling} for smaller NNs (ACAS\_1\_1-student, ACAS\_1\_2-student).
In additon to the smallest and largest MNIST 8x8 NNs considered in previous work (MNIST\_small-top, MNIST\_medium-top), we constructed two
larger MNIST models using student teacher training (MNIST\_large-top, MNIST\_larger-top).
Moreover, we constructed a second version of MNIST\_large-top for $\epsilon$-equivalence verification (MNIST\_large-epsilon).
All NNs were trained using variants of student teacher training and were trained in such a way that they were likely to be top-1 or $\epsilon$-equivalent in some parts of the input space.
More details on the properties of the 9 considered benchmark NNs are available online.\footnote{An overview table of all benchmarks is available at \url{https://github.com/samysweb/nnequiv-experiments/blob/main/benchmarks.md}}
The input space considered for verification is a sensitive choice, as it can have significant and varying impact on the performance of different verification techniques.
For the case of GPE, the algorithm's performance tends to degrade with increasing input space size due to the growth in necessary splits.
Therefore, for each individual benchmark, we decided to look at an input size which was hard to handle for \textsc{NNEquiv-E}.
This has two reasons.
First, it allows us to evaluate the impact of the optimizations presented above in their ability to decrease runtimes.
Secondly, it permits to compare the performance of \textsc{NNEquiv} to the performance of \textsc{MilpEquiv} on instances which are difficult \textit{for our} approach.
The entire experimental setup can be found online.\footnote{On GitHub: \url{https://github.com/samysweb/nnequiv-experiments}}

We used a machine with 4 AMD EPYC 7281 16-Core processors (i.e. 64 cores in total) and a total of 64GBs of RAM.
All experiments were run with a single thread, a memory limit of 512MB\footnote{The memory limit was irrelevant in practice, as no experiment hit this limit.}, and a timeout of 3 hours.
The experiments were run in parallel, up to 24 processes at once.
All times given in the subsequent sections are the median of 3 runs.

\subsection{Comparison of \textsc{NNEquiv} versions}
\label{subsec:experiments-versions}

\Cref{fig:version-runtimes} shows that the proposed optimizations help in reducing the runtime of the algorithm (note that the upper half of the y-axis has a logarithmic scale for improved visibility of the results).
On the one hand, we can observe, that heuristics for $\min_\mathcal{E}$ can, in principle, improve and worsen the result of the approach (as seen with \textsc{NNEquiv-L} and \textsc{NNEquiv-M}).
On the other hand, we see that there is still significant room for improvement through the development of better refinement heuristics -- this optimization would be supplementary to further optimizations which could be developed.

\subsection{Comparison to previous work}
\begin{table}[]
    \centering
    \caption{Runtime comparison (in seconds) for
    \textsc{NNEquiv-L} and \textsc{MilpEquiv}}
    \begin{tabular}{l|c|c|c}
    Benchmark & Property & \textsc{NNEquiv-L} &  \textsc{MilpEquiv} \\
    \hline\hline
    ACAS\_1\_1-retrain &$\epsilon=0.05$&   \textbf{167.45} &       TO \\
    ACAS\_1\_1-student &$\epsilon=0.05$&    \textbf{84.85} &       TO \\
    ACAS\_1\_2-retrain &$\epsilon=0.05$&   \textbf{326.59} &       TO \\
    ACAS\_1\_2-student &$\epsilon=0.05$&   \textbf{109.46} &         320.07 \\
 MNIST\_large-epsilon &$\epsilon=15$&    35.90 &          \textbf{19.97} \\
     MNIST\_small-top &top-1&    14.39 &           \textbf{3.51} \\
    MNIST\_medium-top &top-1&    94.51 &           \textbf{3.85} \\
     MNIST\_large-top &top-1&    \textbf{13.02} &          25.85 \\
    MNIST\_larger-top &top-1&   706.56 &         \textbf{386.04} \\
    \end{tabular}
    \label{tab:nnequiv-milpequiv}
\end{table}
The comparison to \textsc{MilpEquiv} is shown in \Cref{tab:nnequiv-milpequiv}.
\textsc{NNEquiv} outperforms \textsc{MilpEquiv} on the ACAS instances, where \textsc{MilpEquiv} even runs into a timeout for three of the four verification tasks.
In particular, this seems to be the case for larger NNs with low-dimensional inputs.
The superior performance of \textsc{MilpEquiv} for the case of MNIST\_large-epsilon seems to be caused by the LP solver in \textsc{NNEquiv} which is a magnitude slower for solving optimizations tasks for MNIST in comparison to ACAS Xu.
As this cannot be explained by the number of constraints, we suspect this is a problem related to the larger input dimensionality for the MNIST case (64 inputs in comparison to 5 inputs for ACAS Xu).
The ACAS-retrain NNs have the same structure as the original ACAS NN, allowing us to compare \textsc{NNEquiv} to \textsc{ReluDiff}.
While \textsc{ReluDiff} was able to quickly verify equivalence for the truncated NNs, where the mean absolute weight difference was $\approx 9.37*10^{-5}$, it was significantly slower than our approach on the retrain instances, with a mean weight difference of $\approx 0.48$, for $\epsilon \leq 5$ and even timed out for smaller values of $\epsilon$.
This suggests that the applicability of \textsc{ReluDiff} is not only restricted to structurally similar NNs, but that its performance also heavily depends on small weight differences.
Regarding question (E2) we note that our approach is applicable to a broader class of NNs than \textsc{ReluDiff} and solved instances where both other approaches timed out.

\subsection{Influence of $\epsilon$-equivalence tightness}
\begin{figure}
    \centering
    \includegraphics[width=0.9\columnwidth]{figures/epsilon-tightness-alt.png}
    \caption{Solving time in seconds for $\epsilon$-equivalence with varying $\epsilon$: Solving times increase with tighter properties, however \textsc{NNEquiv-L} outperforms \textsc{MilpEquiv}. \textsc{ReluDiff} is outperformed by \textsc{NNEquiv-L} for retrained NNs.}
    \label{fig:epsilon-tightness}
\end{figure}
Concerning question (E3) we evaluated the performance of the approaches as we vary the tightness of $\epsilon$-equivalence for $\epsilon\in\left[0.005,500\right]$.
Note that we did not prove equivalence for $\epsilon=0.005$ for ACAS\_1\_2-student as we found this NN not to be 0.005-equivalent to the original NN.
Intuitively, a proof for a tighter $\epsilon$ bound will require more work as the approach either needs to refine more over-approximations (in the case of \textsc{NNEquiv-L}) or do further branch-and-bound operations (in the case of \textsc{MILPEquiv}).
We can observe this behavior in \Cref{fig:epsilon-tightness} which plots the runtime of \textsc{MILPEquiv} and \textsc{NNEquiv-L} as we tighten the $\epsilon$ bound.
Taking into account the log scales on both axes, we can observe that \textsc{NNEquiv-L} is at least one magnitude faster in proving equivalence for ACAS Xu NNs for $\epsilon\leq0.05$.
In particular, \textsc{MILPEquiv} produces time-outs for 3 of the 4 considered NNs once $\epsilon\leq0.05$.
We therefore suspect that our approach is better at handling very tight $\epsilon$ constraints in large NNs with low dimensional input.
This could potentially be due to the fact that GPE can use additional NN information (layer structure etc.) for its refinement decisions which is not readily available in the branch-and-bound algorithm in the backend of \textsc{MILPEquiv}.
For comparison, we plotted the performance of \textsc{ReluDiff} on NNs for truncated weights and retrained NNs:
As can be seen in \Cref{fig:epsilon-tightness} the approach by Paulsen et al.~\cite{paulsen2020reludiff} behaves similarly with respect to $\epsilon$ tightness.
Additionally, we see that the approach is less efficient for retrained NNs where the equivalence for $\epsilon\leq0.05$ cannot be established.

\subsection{Finding Counterexamples}
Our technique can also be used to find counterexamples, showing that two NNs are not equivalent at a certain point.
This information can be interesting to further train NNs after a failed equivalence proof.
To this end, we compared the capabilities of \textsc{NNEquiv-L} in counterexample finding with the capabilities of \textsc{MilpEquiv}.
To account for possibly easy instances, we looked at a large number of non-equivalent input spaces for each of our benchmark NNs which we know to be equivalent on other parts of the input space.
To generate counterexamples, we randomly sampled input points and selected the points with differing NN outputs.
Using this technique, we produced 100 distinct non-equivalent input points for each of our 9 benchmarks.
These points were used as center for $L_\infty$-balls which represented our input spaces.

We then evaluated \textsc{NNEquiv-L} and \textsc{MilpEquiv} on the same input space radii as in \Cref{subsec:experiments-versions} with the objective of finding counterexamples.
Since counterexample extraction is much faster than equivalence proofs, we set a timeout of 2 minutes.
We found that \textsc{NNEquiv-L} was significantly faster and extracted a counterexample for 890 of the 900 considered benchmarks, while even a version without expensive $\relu$ bounds computation\footnote{While this optimization step improves performance for equivalence proofs, it may degrade performance for counterexample finding.} executed in the initialization (named \textsc{MilpEquiv-B}) was only able to find 315 counterexamples.
Also, the time per
counterexample was significantly lower for \textsc{NNEquiv-L}, making this approach an interesting technique for retraining NNs via counterexamples.
Looking at the behavior of \textsc{MilpEquiv}'s solver backend, it seems that the reason for our superior performance lies in the time needed by \textsc{MilpEquiv} to find an initial feasible solution.
\textsc{MilpEquiv} first needs to resolve the integer based $\relu$ node encodings, which are automatically resolved by \textsc{NNEquiv-L} through the propagation of sets.
\textsc{NNEquiv-L} has the potential to extract polytopes of non-equivalent input space subsets which could allow for even more efficient sampling.
\begin{table}
    \centering
    \caption{Comparison of counterexample finding capabilities of \textsc{NNEquiv-L}, \textsc{MilpEquiv}  and \textsc{MilpEquiv-B} (no $\relu$ node bounds)}
    \begin{tabularx}{\columnwidth}{X|c|c|c}
         & \textsc{NNEquiv-L} & \textsc{MilpEquiv} & \textsc{MilpEquiv-B} \\\hline\hline
        \#Solved & \textbf{890} & 305 & 315\\
        Time~(incl.TO) & \textbf{3,597s} & 75,989s & 72,515s\\
        Time/Solved~(excl.TO) & \textbf{2.69s} & 14.91s & 7.21s\\
    \end{tabularx}
    \label{tab:my_label}
\end{table}

\subsection{Equivalence on Set Data Structures}
For our equivalence verification approach it is necessary to define an equivalence check \texttt{is\_equiv} which verifies whether two set data structures $\Theta_\mathcal{R},\Theta_\mathcal{T}$ satisfy $\epsilon$-equivalence or top-1 equivalence.
First, we present how $\epsilon$-equivalence with Chebyshev norm $\lVert \cdot \rVert_\infty$ can be proven for zonotopes.
Afterwards, we show how Star Sets can be used to prove $\epsilon$-equivalence and top-1 equivalence.

\paragraph{$\epsilon$-Equivalence}
In order to prove $\epsilon$-equivalence with the Chebyshev norm, we need to bound the maximum deviation between the two NN outputs by $\epsilon$.
That is, given the two output zonotopes $\Psi_\mathcal{R} = \langle c_\mathcal{R}, G_\mathcal{R}, P_\mathcal{T} \rangle$ and $\Psi_\mathcal{T} = \langle c_\mathcal{T}, G_\mathcal{T}, P_\mathcal{T} \rangle$ we want to find the maximal deviation:
\begin{align*}
    &\max_{\alpha \in P_\mathcal{T}} \lVert\left(c_\mathcal{R}+G_\mathcal{T} \alpha\right) - \left(c_\mathcal{T} + G_\mathcal{T} \alpha\right)\rVert_\infty\\
    = &\max_{i} \max_{\alpha \in P_\mathcal{T}} \left|\left(c_\mathcal{R}-c_\mathcal{T}\right)_i + \left(G_\mathcal{R}-G_\mathcal{T}\right)_i \alpha\right|.
\end{align*}
As can be seen by the reformulation above, we can find the maximal deviation over the output by solving optimization problems for each dimension of the \textit{differential zonotope}
\[
\partial\Psi = \langle \left(c_\mathcal{R}-c_\mathcal{T}\right), \left(G_\mathcal{R}-G_\mathcal{T}\right), P_\mathcal{T} \rangle).
\]
Recalling that zonotopes can be optimized with a closed form solution, this enables a quick check for the adherence of the desired $\epsilon$-equivalence property.
However, since zonotopes only approximate the output set, one may need to fall back to the use of Star Sets if equivalence cannot be established using zonotopes.
In this case, we can reuse the same formula from above to obtain a \textit{differential star set} $\partial\Theta$ which is then optimized using LP solving.

\paragraph{Top-1 Equivalence}
For top-1 equivalence there are two possible approaches which both rely on propapagated star sets.
We can reuse the \textsc{MilpEquiv} encoding and employ a MILP solver.
Alternatively, we can use a simplex (LP) solver.
In the latter case we split up the output star set $\Theta_\mathcal{T}$:

For every output dimension $1\leq j \leq O$ we generate a polytope $P_j$.
Additional constraints $r_j \geq r_i ~\forall i \neq j$ ensure that output $r_j$ is the maximum among the outputs of $\mathcal{R}$ in $\langle c_\mathcal{R}, G_\mathcal{R}, P_\mathcal{T} \cap P_j\rangle$.
Note that the union of $P_1$ to $P_O$ covers all of $P_\mathcal{T}$.
We then examine the outputs of $\Theta_j=\langle c_\mathcal{T}, G_\mathcal{T}, P_\mathcal{T}\cap P_j\rangle$ for every $1 \leq j \leq O$.
Since $j$ is always the maximum of $\mathcal{R}$ for this part of the output space, we want to ensure that $j$ is also always the maximum of $\mathcal{T}$.
Therefore, we compute the maximal difference between output dimension $j$ and the other dimensions in $\Theta_j$.
If all of these differences are below $0$, we can guarantee top-1 equivalence.
This procedure produces $\mathcal{O}\left(O\right)$ star sets and $\mathcal{O}\left(O^2\right)$ optimization operations in total.

\subsection{Challenges and Limitations of the approach}
While the techniques outlined above permit a straightforward extension of GPE to multiple NNs, and thus allow achieving equivalence verification, the approach comes with a number of pitfalls which should be avoided.
The most obvious is probably the possibility of exponential growth in the number of star sets.
As previously noted, the exact GPE approach based on star sets splits the star sets on ReLU nodes.
Tran et al.~\cite{tran2019star} rightly note that the observed growth usually drastically falls behind the worst case, however the increase in ReLU nodes through the processing of two NNs at once certainly leads to an increase in necessary splits.
This is particularly the case for ReLU nodes which cut off very similar hyperplanes (such as the two ReLU nodes in a NN at the same position with truncated weights).
This can not only double the work, but it may also lead to precision problems with LP solvers which tend to show problematic behavior when encountering a problem which has a very small feasible set\footnote{In one case the solver would return drastically differing maximum values for the same optimization problem depending on the previous requests or would suddenly deem the problem infeasible.}.
To avoid such numerical problems we thus use 64-bit floats by default and always ensure that feasibility is checked at least once by an exact (i.e. rational) LP solver before a branch is declared infeasible.
While this can mitigate most numerical problems,
the approach is weaker than \textsc{ReluDiff}/\textsc{NeuroDiff} for the specific use case of weight truncation for structurally similar NNs (e.g. truncation from 32-bit to 16-bit floats) -- this is best left to the approach presented by Paulsen et al.~\cite{paulsen2020neurodiff}.
Although these initial improvements help in making GPE for equivalence \textit{possible}, this approach is not yet \textit{scalable}.
Hence, we devote the next section to various optimizations.

\subsection{Zonotope Propagation}
As an initial optimization we reused the zonotope propagation technique presented by Bak et al.~\cite{bak2020improved}, which reduces the number of LP optimizations necessary through a zonotope based bounds computation.
We refer to this first version of the algorithm as \textsc{NNEquiv-E} (for exact).
As can be seen in \Cref{fig:version-runtimes} later on, this approach produces a total runtime of 54,390s on our 9 benchmark instances.

\subsection{Zonotope Over-Approximation}
To further optimize the algorithm we can either reduce the time spent per zonotope or we can try to reduce the number of zonotopes which have to be considered.
In order to achieve the second objective, we can over-approximate certain $\relu$ splits through a methodology first presented by Singh et al.~\cite{Singh2018} and later reused by Bak~\cite{10.1007/978-3-030-76384-8_2}:
The idea is to introduce an additional dimension to the zonotope and use it to over-approximate the $\relu$ node by a parallelogram.
Over-approximation errors accumulate across layers (Bak~\cite{10.1007/978-3-030-76384-8_2} refer to this as \textit{error snowball}).
To make the parallelogram as tight as possible and minimize the over-approximation error, we use the bounds computed through LP solving (instead of the looser zonotope bounds) if there were any exact splits beforehand.
In an abstraction-refinement setting, we would start by propagating over-approximating zonotopes through both NNs and then check, whether the equivalence property can be established.
If the property does not hold, we refine one over-approximated $\relu$ node by splitting the zonotope and propagating the split zonotopes further through the NNs.

In \Cref{fig:gpe-splitting} we compare the share of nodes per layer whose bounds contain the value $0$ for an exact approach in comparison to the propagation of an over-approximating zonotope.
Any such node can be considered a split candidate which could be used for refinement.
Each refinement can then help in reducing the over-approximation error and in establishing the desired property.
As can clearly be seen in the plot, the over-approximation approach produces a lot more split candidates than the exact approach.
Not all of the splits candidates encountered for the over-approximation would actually have to be refined in the worst case.
This is, because many of the split candidates are only artifacts of previous over-approximations.
We refer to these split candidates as \textit{ghost splits}.
These \textit{ghost splits} cannot be easily distinguished from actual, necessary splits.
The only guaranteed non-ghost split, is the first split candidate encountered, while all later split candidates might be artifacts of over-approximation.

Thus, the simplest refinement strategy would be to refine only this node.
We refer to this strategy as \textsc{NNEquiv-F} (First), and it reduces the runtime on our benchmark set to 2,489s (c.f. \Cref{fig:version-runtimes}).
However, this approach still leaves
room for improvement as we explain in \Cref{sec:tree-explore}.
\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{figures/split-neurons-acas.png}
    \caption{
    Comparison of exact star set propagation (Exact) and propagation of an over-approximating zonotope (Overapprox) for $\epsilon$-equivalence over ACAS\_1\_1: Net1 corresponds to $\mathcal{R}$ and Net2 corresponds to $\mathcal{T}$ (Overapprox Net2 is hidden behind Overapprox Net1)}
    \label{fig:gpe-splitting}
\end{figure}

\subsection{LP approximation}
With the introduction of over-approximation we encounter an additional problem:
Splitting hyperplanes are no longer dependent on the input variables \textit{only}, but also depend on the dimensions introduced through the over-approximation.
This raises the question how to handle the additional dimension in the propagated star set:
Since equally increasing the dimensionality of the LP problem leads to increased solver runtimes,
we instead opted to over-approximate the LP problem.
Classically, for an $m$ dimensional zonotope with initial input dimensionality $I$ we observe a hyperplane cut of the following form:
\begin{align}
\sum_{i=1}^I g_i \alpha_i + \sum_{i=I+1}^m g_i \alpha_i \leq c
\label{eq:exact-lp}
\end{align}
We can now over-approximate this inequality by computing $\mu = \min_{\alpha} \sum_{i=I+1}^m g_i \alpha_i$ through zonotope optimization and constraining the LP problem with the following inequality:
\begin{align}
    \sum_{i=1}^I g_i \alpha_i \leq c - \mu
    \label{eq:approx-lp}
\end{align}
Since any solution for \Cref{eq:exact-lp} implies a solution for \Cref{eq:approx-lp} the second inequality is an over-approximation and can be used to reduce the number of dimensions the LP solver has to handle despite the
over-approximation of the zonotope
Note that we need to take this over-approximation into account for minimization/maximization tasks.
Since the LP solver only optimizes the first $I$ dimensions, we need to add the optimization result of the over-approximating zonotope for the remaining dimensions.
We refer to this version as \textsc{NNEquiv-A} (for approximate LP).
\Cref{fig:version-runtimes} shows that this approach
reduces the runtime to 1,631s.

\subsection{Feed-Forward NNs}
NNs consist of interconnected units called neurons.
A neuron $j$ computes a non-linear function of its input values $x_1, \dots, x_n$
according to $y_j = \sigma(\sum_{i=1}^{n} w_{ij} x_i + b_j)$
where $\sigma$ is called the activation function and $w_{ij}$ are the weights.
$b_j$ is commonly referred to as bias of neuron $j$.
In this paper, we focus on the \emph{rectified linear unit} activation function,
$\relu(x) = \max(0, x)$, which is one of the most commonly used activation functions in modern NNs~\cite{Goodfellow-et-al-2016}.
Outputs of neurons are connected as input to other neurons, resulting in a directed graph.
In this paper, we focus on feed-forward NNs, where the underlying graph is acyclic.
Neurons are organized in layers, where neurons in layer $l$ take inputs only from the directly preceding layer $l-1$.
The first layer---called input layer---is just a place holder for the inputs to be fed into the NN, the subsequent layers are called hidden layers, while the last layer---the output layer---holds the function value computed by the NN.
We refer to the input space dimension as $I\in\mathbb{N}$ and to the output dimension as $O\in\mathbb{N}$.

\subsection{Verifiying NN Equivalence}
There has recently been a line of work which proposes various compression techniques for NNs (for a full review see Cheng et al.~\cite{cheng2020survey}).
While such techniques have been shown to be useful in practice, most lack a formal proof of correctness and only  rely on empirical evidence.
The usage of such techniques thus raises the question of how to prove that two NNs $\mathcal{R}$ (reference) and $\mathcal{T}$ (test) and their corresponding mathematical functions $g_{\mathcal{R}}:\mathbb{R}^I\to \mathbb{R}^O$, $g_{\mathcal{T}}:\mathbb{R}^I\to \mathbb{R}^O$ are equivalent, i.e. that they produce the same results.

Kleine Büning et al.~\cite{Buning2020} examined three possible definitions of equivalence
on a given subset of inputs $\mathcal{I}\subseteq \mathbb{R}^I$,
two of which we review here
\begin{definition}[$\epsilon$-Equivalence \cite{Buning2020} / Differential Equivalence~\cite{paulsen2020reludiff}]
\label{def:epsilon-equiv}
Two NNs $\mathcal{R}$ and $\mathcal{T}$ are $\epsilon$-equivalent with respect to a norm $\lVert \cdot \rVert$, if $\left\lVert g_{\mathcal{R}}\left(x\right)-g_{\mathcal{T}}\left(x\right)\right\rVert < \epsilon$ for all $x \in \mathcal{I}$
\end{definition}
\begin{definition}[Top-1-Equivalence \cite{Buning2020}]
\label{def:one-hot-equiv}
Two NNs $\mathcal{R}$ and $\mathcal{T}$ are top-1-equivalent, if $\argmax_i r_i = \argmax_j t_j$ where $r=g_{\mathcal{R}}\left(x\right)$ and $t=g_{\mathcal{T}}\left(x\right)$ for all $x \in \mathcal{I}$
\end{definition}
While Definition \ref{def:epsilon-equiv} is particularly suitable for regression tasks as it can show a very strong form of equivalence, the latter definition is more relaxed and especially useful for classification tasks.

Paulsen et al.~\cite{paulsen2020neurodiff} and Kleine Büning et al.~\cite{Buning2020} proposed two fundamentally different approaches to NN equivalence verification.
While Paulsen et al.~\cite{paulsen2020neurodiff} proposed a technique called \textsc{ReluDiff}/\textsc{NeuroDiff} which uses (symbolic) interval propagation on NNs with similar weight configurations (e.g. produced through float truncation)
to prove $\epsilon$-equivalence, the \textsc{MilpEquiv} technique \cite{Buning2020} is based on mixed integer linear programming (MILP) and encodes (potentially) structurally different NNs together with the desired equivalence property into an optimization problem.
\textsc{ReluDiff}/\textsc{NeuroDiff} was shown to be efficient for cases where a NN's weights had been truncated.
However, this approach does not work at all for structurally differing NNs and suffers in performance when weight differences are larger.
For such NNs \textsc{MilpEquiv} was shown to work well for a number of small instances.
(e.g. two 64 pixel input MNIST NNs with a total of 84 $\relu$ nodes).
In particular, \textsc{MilpEquiv} is able to provide a maximal radius around a datapoint for which the property still holds.

\subsection{Geometric Path Enumeration}
\label{subsec:related-gpe}
GPE is a methodology originally proposed by Tran et al.~\cite{tran2019star} for verifying safety properties in NNs.
Given a NN $\mathcal{N}$, a set of input instances $\mathcal{I}\subseteq\mathbb{R}^I$ and an unsafe output specification $\mathcal{U}\subseteq\mathbb{R}^O$ defined as a set of linear constraints,
safety verification is concerned with the question whether there exist any instances $i\in\mathcal{I}$ such that $g_{\mathcal{N}}\left(i\right) \in \mathcal{U}$.

Instead of pushing single data points through the NN and checking whether they satisfy the required safety property, GPE feeds an entire set into the NN and then evaluates whether any parts of the output sets lie inside $\mathcal{U}$.
The sets are represented through generalized star sets or zonotopes which we define below:

\begin{definition}[Generalized Star Set \cite{tran2019star}]
A generalized star set $\Theta$ is a tuple $\left\langle c, G, P\right\rangle$ where $c\in\mathbb{R}^n$ is the center, $G=\left(g_1\cdots g_m\right)\in\mathbb{R}^{n \times m}$ is the generator matrix, and $P\subseteq\mathbb{R}^m$ is a set defined through a conjunction of linear constraints (i.e. a polytope).
The set represented by $\Theta$ is then defined as:
\[
\left\llbracket \Theta \right\rrbracket = \left\{ x \in \mathbb{R}^n \;\middle\vert\; \exists \alpha \in P \,:\, x = c + G\alpha \right\}
\]
\label{def:gpe:star-set}
\end{definition}
Alternatively, it is possible to use zonotopes which further restrict the type of predicate allowed:
\begin{definition}[Zonotopes \cite{bak2020improved}]
A zonotope $\Psi$ is a generalized star set $\langle c, G, P \rangle$ with the further restriction that $P$ may only be defined through interval constraints (i.e. $P$ only enforces a lower and upper bound for each dimension).
\end{definition}

We refer to zonotopes and generalized star sets as \emph{set data structures}.
Initially, the GPE algorithm converts the provided input space into either of these data structures and then propagates the sets through the NN.
The transformation deployed through the dense layer of a NN can be exactly represented by zonotopes and generalized star sets through application of the weight matrices to $G$ and $c$.
$\relu$ nodes require a different type of transformation since they are only piece-wise linear functions:
To this end, the data structures can either be split by introduction of an additional hyperplane to the linear constraint predicate (\textit{exact GPE}) or the $\relu$ function can be over-approximated~\cite{10.1007/978-3-030-76384-8_2,Singh2018} (\textit{approximate GPE}).
Note that optimization for generalized star sets is much more expensive than optimization for zonotopes, which can be computed with a closed form solution.