Probabilistic verification of neural networks is concerned with formally analysing the output distribution of a neural network under a probability distribution of the inputs.
Examples of probabilistic verification include verifying the demographic parity fairness notion or quantifying the safety of a neural network.
We present a new algorithm for the probabilistic verification of neural networks based on an algorithm for computing and iteratively refining lower and upper bounds on probabilities over the outputs of a neural network.
By applying state-of-the-art bound propagation and branch and bound techniques from non-probabilistic neural network verification, our algorithm significantly outpaces existing probabilistic verification algorithms, reducing solving times for various benchmarks from the literature from tens of minutes to tens of seconds.
Furthermore, our algorithm compares favourably even to dedicated algorithms for restricted subsets of probabilistic verification.
We complement our empirical evaluation with a theoretical analysis, proving that our algorithm is sound and, under mildly restrictive conditions, also complete when using a suitable set of heuristics.

§ INTRODUCTION

As deep learning spreads through society, it becomes increasingly important to ensure the reliability of artificial neural networks, including aspects of fairness and safety.

However, manually introspecting neural networks is infeasible due to their opaque nature, and empirical assessments of neural networks are challenged by neural networks being fragile with respect to various types of input perturbations <cit.>.
In contrast, neural network verification analyses neural networks with mathematical rigour, facilitating the faithful auditing of neural networks.

In this paper, we consider probabilistic verification of neural networks, which is concerned with proving statements about the output distribution of a neural network.

An example of probabilistic verification is proving that a neural network \(\NN\) making a binary decision affecting a person (for example, hire/do not hire, credit approved/denied) satisfies the demographic parity fairness notion <cit.> under a probability distribution \(\prob[\inp]\) of the network inputs \(\inp\) representing the person

\begin{equation}
\frac{%
\prob[\vec{x}](\NN(\vec{x}) = \texttt{yes} \mid \vec{x}\ \text{is disadvantaged})
\prob[\vec{x}](\NN(\vec{x}) = \texttt{yes} \mid \vec{x}\ \text{is advantaged})
} \geq \gamma,\label{eqn:demographic-parity-intro}
\end{equation}

where \(\gamma \in [0, 1]\) with \(\gamma = 0.8\) being a common choice <cit.>.
A closely related problem to probabilistic verification is computing bounds on probabilities over a neural network.
An example of this is quantifying the safety of a neural network by bounding

\begin{equation}
\prob[\vec{x}](\NN(\vec{x})\ \text{is unsafe}).\label{eqn:safety-intro}
\end{equation}

In this paper, we introduce a novel algorithm for computing bounds on probabilities such as <ref> using a branch and bound framework <cit.>.
These bounds then allow us to verify probabilistic statements like <ref> using ideas from the probabilistic verification algorithm  <cit.>.

More concretely, we recombine neural network verification using branch and bound <cit.>, linear relaxations of neural networks <cit.>, and massively parallel branch and bound <cit.> with ideas from  <cit.> for probabilistic neural network verification, to obtain the * algorithm, a fast and generally applicable probabilistic verification algorithm for neural networks.
Our theoretical analysis of  shows that  is sound and, under mildly restrictive conditions, complete when using suitable branching and splitting heuristics.

Our experimental evaluation reveals that  significantly outpaces the probabilistic verification algorithms  <cit.> and  <cit.>.
In particular, we solve benchmark instances that  can not solve within 15 minutes in less than 10 seconds and solve the ACAS Xu <cit.> probabilistic robustness case study of ConverseFilieriGopinathEtAl2020 in a mean runtime of 18 seconds, compared to 33 minutes for .

Applying  to #DNN verification <cit.>, a subset of probabilistic verification, reveals that  also compares favourably to the  algorithm <cit.> specialised to #DNN verification and even the  algorithm <cit.> that relaxes #DNN verification to computing a confidence interval on the solution.
In contrast to  and similar approaches <cit.>,  computes lower and upper bounds on probabilities like <ref> that are guaranteed to hold with absolute certainty.
Such bounds are preferable to confidence intervals in high-risk machine-learning applications.

To test the limits of  and fuel further research in probabilistic verification, we introduce a significantly more challenging probabilistic verification benchmark:
MiniACSIncome is a benchmark based on the ACSIncome dataset <cit.> and is concerned with verifying the demographic parity of neural networks for datasets of increasing input dimensionality.

In summary, our contributions are

* the  algorithm for the probabilistic verification of neural networks,

* a theoretical analysis of ,

* a thorough experimental comparison of  with existing probabilistic verifiers for neural networks and tools dedicated to restricted subsets of probabilistic verification, and

* MiniACSIncome: a new, challenging probabilistic verification benchmark.

§ RELATED WORK

Non-probabilistic neural network verification is concerned with proving that the outputs of a neural network satisfy some condition for all inputs in an input set.
Approaches for non-probablistic neural network verification include Satisfiability Modulo Theories (SMT) solving <cit.>, Mixed Integer Linear Programming (MILP) <cit.>, and Reachability Analysis <cit.>.
Many of these approaches can be understood as branch and bound algorithms <cit.>.
Branch and bound <cit.> also powers the \(\alpha,\!\beta\)-CROWN <cit.>, MN-BaB <cit.> and VeriNet <cit.> verifiers that lead the table in recent international neural network verifier competitions <cit.>.

A critical component of a branch and bound verification algorithm is computing bounds on the output of a neural network.
Approaches for bounding neural network outputs include interval arithmetic <cit.>, dual approaches <cit.>, linear bound propagation techniques <cit.>, multi-neuron linear relaxations <cit.>, and further optimisation-based approaches <cit.>.

Probabilistic verification algorithms can be divided into sound algorithms that provide valid proofs and probably sound algorithms that provide valid proofs with a certain predefined probability.
<cit.> is a fairness verification algorithm, a subset of probabilistic verification that studies problems such as <ref>.
uses SMT solving for splitting the input space into disjoint hyperrectangles, integrating over which yields bounds on a target probability.
ConverseFilieriGopinathEtAl2020 () and BorcaTasciucGuoBakEtAl2023 divide the input space into disjoint polytopes using concolic execution and reachable set verification, respectively, to perform probabilistic verification.
MorettinPasseriniSebastiani2024 use weighted model integration <cit.> to obtain a general probabilistic verification algorithm but neither provide code nor report runtimes.
A restricted subset of probabilistic verification is #DNN verification <cit.>, corresponding to probabilistic verification under uniformly distributed inputs.
The #DNN verifier  <cit.> uses a similar massively parallel branch and bound approach as our algorithm but does not perform general probabilistic verification.

Probably sound verification algorithms <cit.> obtain efficiency at the cost of potentially unsound results.
<cit.> is a probably sound #DNN verification algorithm.
ConverseFilieriGopinathEtAl2020,BastaniZhangSolarLezama2019,MarzariRoncolatoFarinelli2023 compare sound and probably sound approaches for fairness verification, general probabilistic verification, and #DNN verification, respectively.
We study sound probabilistic verification since certainly sound results are preferable in critical applications, such as education <cit.>, medical applications, or autonomous driving and flight.

Besides probabilistic fairness notions, such as demographic parity <cit.>, several approaches <cit.> verify dependency fairness <cit.>, an individual fairness notion <cit.> that states that persons that only differ by their protected attribute (for example, gender or race) need to be assigned to the same class.
However, dependency fairness can be satisfied trivially by withholding the protected attribute from the classifier.
Since withholding the protected attribute is insufficient <cit.> or even harmful <cit.> for fairness, dependency fairness is an insufficient fairness notion.

§ PRELIMINARIES AND PROBLEM STATEMENT

Throughout this paper, we are concerned with computing (provable) lower and upper bounds on various functions.

For \(f: \Reals^n \to \Reals^m\), we call \(\ell, u \in \Reals^m\) a lower, respectively, upper bound on \(f\) for \(\mathcal{X}' \subseteq \Reals^n\) if \(\ell \leq f(\vec{x}) \leq u, \forall \vec{x} \in \mathcal{X}'\).

Neural Networks.
In particular, we are concerned with computing bounds on functions involving neural networks \(\NN: \mathcal{X} \to \Reals^m\), where \(\mathcal{X} \subset \Reals^n\) is the input space of the neural network.
A neural network is a composition of linear functions and a predefined set of non-linear functions, such as ReLU, Tanh, and max pooling.
We refrain from further defining neural networks but refer the interested reader to the auto_LiRPA library <cit.> that practically defines the class of neural networks to which our algorithm can be applied.
For the scope of this paper, we only consider fully-connected feed-forward neural networks.

Throughout this paper, we assume that the input space \(\mathcal{X}\) is a bounded hyperrectangle.
Bounded input spaces are common in machine learning.
Examples include the space of normalised images \({[0, 1]}^n\) in computer vision and tabular input spaces of practically bounded variables, such as age, working hours per week, income, etc.

Notation and Terminology.

We use bold letters (\(\vec{x}\), \(\vec{y}\)) for vectors and calligraphic letters (\(\mathcal{X}\)) for sets.
Let \(\underline{\vec{x}}, \overline{\vec{x}} \in \Reals^n\) with \(\underline{\vec{x}} \leq \overline{\vec{x}}\).
We use \([\underline{\vec{x}}, \overline{\vec{x}}] = \{\vec{x} \in \Reals^n \mid \underline{\vec{x}} \leq \vec{x} \leq \overline{\vec{x}}\}\) to denote the hyperrectangle with the minimal element \(\underline{\vec{x}}\) and the maximal element \(\overline{\vec{x}}\).
When we speak of bounds on a certain quantity in this paper, we refer to a pair of a lower and an upper bound.
Throughout this paper, \(\ell, \ka\) denote lower bounds, and \(u\) denotes an upper bound.
When convenient, we also use \(\underline{\vec{y}} \leq \vec{y} \leq \overline{\vec{y}}\) to denote bounds, especially on vectors.

We assume that all random objects are defined on the same abstract probability space $(\Omega, \mathcal{F}, \prob)$ and that all continuous random variables admit a probability density function.

§.§ Probabilistic Verification of Neural Networks

Let \(\mathcal{X} \subset \Reals^n\) be a bounded hyperrectangle and let \(v \in \Nats\).
In this paper, we are concerned with proving or disproving whether \(\NN: \mathcal{X} \to \Reals^m\) is feasible for the probabilistic verification problem

\begin{equation}
\PVP: \left\{\,
\begin{array}{l}
\fSAT(p_1, \ldots, p_v) \geq 0,\\
p_i = \prob[\inp[i]](
\fSAT*[i](\inp[i], \NN(\inp[i])) \geq 0
) \quad \forall i \in \{1, \ldots, v\},
\end{array}
\right.\label{eqn:pvp}
\end{equation}

where \(\inp[i]\), \(i \in \{1, \ldots, v\}\), is a $\mathcal{X}$-valued random variable with distribution \(\prob[\inp[i]]\) and \(\fSAT: \Reals^v \to \Reals\), \(\fSAT*[i]: \Reals^{n} \times \Reals^{m} \to \Reals\), \(i \in \{1, \ldots, v\}\) are satisfaction functions.
We only consider satisfaction functions that are compositions of linear functions, multiplication, division, and monotone functions, such as ReLU, Sigmoid, \(\min\) and \(\max\).
Throughout this paper, we assume all probabilistic verification problems to be well-defined.

We express the demographic parity fairness notion from <ref> as a probabilistic verification problem.

Let \(\mathcal{X} \subset \Reals^n\) be an input space containing a categorical protected attribute, such as gender, race, or disability status that is one-hot encoded at the indices \(A \subset \{1, \ldots, n\}\).
We assume a single historically advantaged category encoded at the index \(a \in A\).

Consider a neural network \(\NN: \Reals^n \to \Reals^2\) that acts as a binary classifier making a decision affecting a person, such as hiring or credit approval.
The neural network produces a score for each class and assigns the class with the higher score to an input.

We express <ref> as a probabilistic verification problem, that is,

\begin{equation*}
\frac{%
\prob[\vec{x}](\NN(\inp) = \texttt{yes} \mid \inp\ \text{is disadvantaged})
\prob[\vec{x}](\NN(\inp) = \texttt{yes} \mid \inp\ \text{is advantaged})
\geq \gamma
\quad \Longleftrightarrow \quad
\NN\ \text{is feasible in \cref{eqn:pvp}},
\end{equation*}

where, \(\fSAT(p_1,p_2,p_3,p_4) = (p_1 p_4)/(p_2 p_3) - \gamma\), \(\fSAT*[1](\inp, \NN(\inp)) = \min({\NN(\inp)}_1 - {\NN(\inp)}_2, -\inp_a)\), \(\fSAT*[2](\inp, \NN(\inp))=-\inp_a\), \(\fSAT*[3](\inp, \NN(\inp)) = \min({\NN(\inp)}_1 - {\NN(\inp)}_2, \inp_a - 1)\), and \(\fSAT*[4](\inp, \NN(\inp)) = \inp_a -1\).
<Ref> contains a detailed derivation of this equivalence.

If all inputs \(\inp[i]\) are uniformly distributed, probabilistic verification corresponds to #DNN verification <cit.>.
As MarzariCorsiCicaleseEtAl2023 prove, #DNN verification is #P complete, implying that probabilistic verification is #P hard.
Since #P is at least as hard as NP, we can not expect to obtain efficient algorithms, but which probabilistic verification problems are practically solvable remains undetermined.

§.§ Non-Probabilistic Neural Network Verification

The goal of non-probabilistic neural network verification is to prove or disprove whether a neural network \(\NN: \mathcal{X} \to \Reals^m\) is feasible for

\begin{equation}
\VP: \;%\Big\{
\parg{} \quad \forall \vec{x} \in \mathcal{X}',\label{eqn:nn-verif}
\end{equation}

where \(\mathcal{X}' \subseteq \mathcal{X}\) is a hyperrectangle, and \(\fSAT*:\Reals^{n} \times \Reals^{m} \to \Reals\) is a satisfaction function that indicates whether the output of \(\NN\) is desirable (\(\fSAT*(\cdot, \NN(\cdot)) \geq 0\)) or undesirable (\(\fSAT*(\cdot, \NN(\cdot)) < 0\)).
In (non-probabilistic) neural network verification, \(\fSAT*\) can generally be considered a part of \(\NN\) <cit.>.

Neural network verifiers are algorithms for proving or disproving <ref>.
Two desirable properties of neural network verifiers are soundness and completeness.

A verification algorithm is sound if it only produces genuine counterexamples and valid proofs for <ref>.
It is complete if it produces a counterexample or proof for <ref> for any neural network in a finite amount of time.

Analogous notions of soundness and completeness also apply to probabilistic verifiers.
Generally unsound but probably sound probabilistic verification approaches are discussed in <ref>.
We introduce incomplete and complete neural network verification through two examples: interval arithmetic and branch and bound.

§.§.§ Interval Arithmetic

Interval arithmetic <cit.> is a bound propagation technique that derives bounds on the output of a neural network from bounds on the network input.

Assume \(\underline{\inp} \leq \inp \leq \overline{\inp}\) are bounds on the network input \(\inp\) and we apply interval arithmetic to compute \(\ell \leq \fSAT*(\inp, \NN(\inp)) \leq u\), \(\forall \inp \in \HR{\inp}\).
If the lower bound is large enough or the upper bound small enough, we can prove or disprove <ref> using

\begin{equation*}
\ell \geq 0 \Longrightarrow \fSAT*(\inp, \NN(\inp)) \geq 0,
\qquad\text{respectively,}\qquad
u < 0 \Longrightarrow \fSAT*(\inp, \NN(\inp)) < 0.
\end{equation*}

However, if the bounds are inconclusive, that is \(\ell < 0 \leq u\), we can neither prove nor disprove <ref>.
The possibility of inconclusive results makes interval arithmetic an incomplete technique for neural network verification.
In particular, interval arithmetic is relatively imprecise, meaning that bounds derived using interval arithmetic are frequently inconclusive <cit.>.
More precise approaches include DeepPoly <cit.> and  <cit.> that follow the same overall scheme but propagate linear bounds through a neural network.
We primarily employ interval arithmetic to compute bounds on the \(\fSAT\) function in <ref> given bounds on \(p_1, \ldots, p_v\).

Let \(f: \Reals^n \to \Reals\) be a function that we want to bound for inputs \(\inp \in \HR{\inp}\).
Interval arithmetic and other bound propagation techniques rely on \(f = f^{(K)} \circ \cdots \circ f^{(1)}\) being a composition of more fundamental functions \(f^{(k)}: \Reals^{n_k} \to \Reals^{n_{k+1}}\) for which we can already compute \(\ell^{(k)} \leq f^{(k)}(\vec{z}) \leq u^{(k)}\) given \(\underline{\vec{z}} \leq \vec{z} \leq \overline{\vec{z}}\).
Examples of such functions include monotone non-decreasing functions, for which \(f^{(k)}(\underline{\vec{z}}) \leq f^{(k)}(\vec{z}) \leq f^{(k)}(\overline{\vec{z}})\), or monotone non-increasing functions, for which \(f^{(k)}(\overline{\vec{z}}) \leq f^{(k)}(\vec{z}) \leq f^{(k)}(\underline{\vec{z}})\).
These two bounding rules already allow us to bound addition, subtraction, \(\min\), \(\max\), and many popular neural network activation functions like ReLU, Sigmoid, and Tanh.
<Ref> contains bounding rules for linear functions, multiplication, and division.
<Ref> describes how interval arithmetic computes bounds on \(f\) using the bounding rules for \(f^{(K)} \ldots, f^{(1)}\).

Interval Arithmetic
Function \(f = f^{(K)} \circ \ldots \circ f^{(1)}\), Input Bounds \(\underline{\inp}, \overline{\inp}\)
\(\underline{\vec{z}}^{(0)} \gets \underline{\inp}\), \(\overline{\vec{z}}^{(0)} \gets \overline{\inp}\)
\(k \in \{1, \ldots, K\}\)
\((\underline{\vec{z}}^{(k)}, \overline{\vec{z}}^{(k)}) \gets \BoundingRule(f^{(k)}, \underline{\vec{z}}^{(k-1)}, \overline{\vec{z}}^{(k-1)})\)

\(\underline{\vec{z}}^{(K)}, \overline{\vec{z}}^{(K)}\)

§.§.§ Branch and Bound

Irrespective of their relative precision, bound propagation approaches such as interval arithmetic, DeepPoly <cit.>, and * are incomplete according to <ref>.
To obtain a complete verifier, bound propagation can be combined with branching to gain completeness.
This algorithmic framework is called branch and bound <cit.>.
In branch and bound, when the computed bounds are inconclusive (\(\ell < 0 \leq u\)), the search space is split (branching).
The idea is that splitting improves the precision of the bounds for each part of the split (each branch).

<Ref> contains an abstract branch and bound algorithm for neural network verification.
This algorithm has three subprocedures that we need to instantiate for running <ref>: , , and .
One instantiation that yields a complete algorithm is selecting branches in a FIFO order, computing bounds using interval arithmetic and splitting the search space by bisecting the input space <cit.>.
XuZhangWangEtAl2021 present an improved complete branch and bound algorithm that selects branches in batches to utilise massively parallel hardware such as GPUs, computes bounds using a combined bound propagation and optimisation procedure, and splits the search space by splitting on hidden ReLU nodes in the network <cit.>.
We refer to ZhangWangXuEtAl2022 for a state-of-the-art neural network verifier <cit.> based on branch and bound.
In the next section, we introduce a branch and bound algorithm for the probabilistic verification of neural networks.

Branch and Bound for Non-Probabilistic Neural Network Verification
Verification Problem \(\fSAT*(\inp, \NN(\inp)) \geq 0, \forall \vec{x} \in \mathcal{X}'\) where \(\mathcal{X}' = \HR{\inp}\)

\(\mathrm{branches} \gets \{\mathcal{X}'\}\)
\(\mathrm{branches} \neq \emptyset\)
\(\HR[']{\vec{x}} \gets \Select(\mathrm{branches})\)*Select one branch
\((\underline{y}, \overline{y}) \gets \ComputeBounds(\fSAT*(\cdot, \NN(\cdot)), \underline{\inp}, \overline{\inp})\)*E.g. interval arithmetic
(*[f]Every \(\vec{x} \in \HR[']{\vec{x}}\) is a counterexample)\(\overline{y} < 0\)

(*[f]Prune branch)\(\underline{y} \geq 0\)
\(\mathrm{branches} \gets (\mathrm{branches} \setminus \{\HR[']{\vec{x}}\})\)

\(\mathrm{new} \gets \Split(\HR[']{\vec{x}})\)*E.g. by bisection
\(\mathrm{branches} \gets (\mathrm{branches} \setminus \{\HR[']{\vec{x}}\}) \cup \mathrm{new}\)

*All branches pruned

§ ALGORITHM

This section introduces *, our algorithm for probabilistic verification of neural networks as defined in <ref>.

For conciseness, we use \(p_i\) for \(i \in \{1, \ldots, v\}\) as in <ref> and drop the index and superscripts \({}^{(i)}\) whenever we only consider a single probability \(p\).

uses the same overall approach as  <cit.>: Iteratively refine bounds on each \(p_i\) until interval arithmetic allows us to prove or disprove \(\fSAT(p_1, \ldots, p_v) \geq 0\) (see <ref>).
We also follow  in splitting the input space into hyperrectangles since this allows for computing probabilities efficiently.
However, while  uses expensive SMT solving for refining the input space, we use a branch and bound algorithm utilising computationally inexpensive input splitting and bound propagation techniques from non-probabilistic neural network verification for refining the input splitting.
<Ref> illustrates our approach for computing and refining bounds on a probability \(p\).
The following sections first introduce the algorithm that solves <ref> given bounds on \(p_1, \ldots, p_v\) before introducing the algorithm for computing the bounds on each \(p_i\).

side by side/.style 2 args=
line width=8pt,

class A/.style=
class B/.style=
class unknown/.style=

axis equal image,
scale only axis,
xlabel=$\vec{x}_1$, ylabel=$\vec{x}_2$,
xmin=0, xmax=10,
ymin=0, ymax=10,
xtick=, ytick=,

(D0)  at (-1,7.0);
(D1)  at (-0.1,7.45);
(D2)  at (1.6,8.2);
(D3)  at (4.0,8.5);
(D4)  at (4.6,8.1);
(D5)  at (4.75,7.3);
(D6)  at (4.65,6.3);
(D7)  at (4.4,4.7);
(D8)  at (4.5,3.8);
(D9)  at (4.7,2.6);
(D10) at (5.1,2.0);
(D11) at (5.5,1.3);
(D12) at (5.8,-0.2);
(D13) at (6.0,-1);

(S1_0) at (5,0);
(S1_1) at (5,10);

(S2_0) at (0,5);
(S2_1) at (5,5);

(S3_0) at (2.5,0);
(S3_1) at (2.5,5);

(S4_0) at (5,5);
(S4_1) at (10,5);

[fill=Class1] (0,0) – (0,5) – (2.5,5) – (2.5,0) – cycle;
[fill=Class2] (5,5) – (5,10) – (10,10) – (10,5) – cycle;
[fill=black!10!white] (0,5) – (0,10) – (5,10) – (5,5) – cycle;

[fill=black!10!white] (2.5,0) – (2.5,5) – (5,5) – (5,0) – cycle;
[fill=black!10!white] (5,0) – (5,5) – (10,5) – (10,0) – cycle;

[gray,thin,step=1.25] (0,0) grid (10,10);

[side by side=Class2Class1] (D0) – (D1) – (D2) – (D3) – (D4) – (D5) – (D6) – (D7) – (D8) – (D9) – (D10) – (D11) – (D12) – (D13);
[black,thick] (D0) – (D1) – (D2) – (D3) – (D4) – (D5) – (D6) – (D7) – (D8) – (D9) – (D10) – (D11) – (D12) – (D13);

[black,very thick]
(S1_0) – (S1_1)
(S2_0) – (S2_1)
(S3_0) – (S3_1)
(S4_0) – (S4_1)

[class A] at (1.25,2.5) V;
[class B] at (7.5,7.5) S;

[pin=[pin edge=black,semithick,Stealth-,font=,align=center,xshift=-1.1cm,yshift=-2.05cm]Satisfaction
Boundary] at ($(D4) + (.09cm,.09cm)$) ;

After four Splits.

axis equal image,
scale only axis,
xlabel=$\vec{x}_1$, ylabel=$\vec{x}_2$,
xmin=0, xmax=10,
ymin=0, ymax=10,
xtick=, ytick=,

(D0)  at (-1,7.0);
(D1)  at (-0.1,7.45);
(D2)  at (1.6,8.2);
(D3)  at (4.0,8.5);
(D4)  at (4.6,8.1);
(D5)  at (4.75,7.3);
(D6)  at (4.65,6.3);
(D7)  at (4.4,4.7);
(D8)  at (4.5,3.8);
(D9)  at (4.7,2.6);
(D10) at (5.1,2.0);
(D11) at (5.5,1.3);
(D12) at (5.8,-0.2);
(D13) at (6.0,-1);

(S1_0) at (5,0);
(S1_1) at (5,10);

(S2_0) at (0,5);
(S2_1) at (5,5);

(S3_0) at (2.5,0);
(S3_1) at (2.5,5);

(S4_0) at (5,5);
(S4_1) at (10,5);

(S6_0) at (7.5,0);
(S6_1) at (7.5,5);

(S7_0) at (5,2.5);
(S7_1) at (10,2.5);

(S8_0) at (5,2.5);
(S8_1) at (10,2.5);

(S9_0) at (6.25,0);
(S9_1) at (6.25,2.5);

(S11_0) at (2.5,2.5);
(S11_1) at (5,2.5);

(S12_0) at (3.75,2.5);
(S12_1) at (3.75,5);

(S13_0) at (2.5,5);
(S13_1) at (2.5,10);

(S14_0) at (0,7.5);
(S14_1) at (2.5,7.5);

(S15_0) at (2.5,7.5);
(S15_1) at (5,7.5);

(S16_0) at (3.75,0);
(S16_1) at (3.75,2.5);

(S17_0) at (3.75,5);
(S17_1) at (3.75,7.5);

(S18_0) at (2.5,8.75);
(S18_1) at (5,8.75);

(S19_0) at (0,8.75);
(S19_1) at (2.5,8.75);

[fill=Class1] (0,0) – (0,5) – (2.5,5) – (2.5,0) – cycle;
[fill=Class1] (2.5,2.5) – (2.5,5) – (3.75,5) – (3.75,2.5) – cycle;
[fill=Class1] (2.5,0) – (2.5,2.5) – (3.75,2.5) – (3.75,0) – cycle;
[fill=black!10!white] (3.75,2.5) – (3.75,5) – (5,5) – (5,2.5) – cycle;
[fill=black!10!white] (3.75,0) – (3.75,2.5) – (5,2.5) – (5,0) – cycle;

[fill=Class1] (0,5) – (0,7.5) – (2.5,7.5) – (2.5,5) – cycle;
[fill=Class1] (2.5,5) – (2.5,7.5) – (3.75,7.5) – (3.75,5) – cycle;
[fill=Class2] (0,8.75) – (0,10) – (2.5,10) – (2.5,8.75) – cycle;
[fill=Class2] (2.5,8.75) – (2.5,10) – (5,10) – (5,8.75) – cycle;
[fill=black!10!white] (0,7.5) – (0,8.75) – (2.5,8.75) – (2.5,7.5) – cycle;
[fill=black!10!white] (2.5,7.5) – (2.5,8.75) – (5,8.75) – (5,7.5) – cycle;
[fill=black!10!white] (3.75,5) – (3.75,7.5) – (5,7.5) – (5,5) – cycle;

[fill=Class2] (5,5) – (5,10) – (10,10) – (10,5) – cycle;

[fill=Class2] (5,2.5) – (5,5) – (7.5,5) – (7.5,2.5) – cycle;
[fill=Class2] (7.5,2.5) – (7.5,5) – (10,5) – (10,2.5) – cycle;
[fill=Class2] (7.5,0) – (7.5,2.5) – (10,2.5) – (10,0) – cycle;
[fill=Class2] (6.25,0) – (6.25,2.5) – (7.5,2.5) – (7.5,0) – cycle;
[fill=black!10!white] (5,0) – (5,2.5) – (6.25,2.5) – (6.25,0) – cycle;

[gray,thin,step=1.25] (0,0) grid (10,10);

[side by side=Class2Class1] (D0) – (D1) – (D2) – (D3) – (D4) – (D5) – (D6) – (D7) – (D8) – (D9) – (D10) – (D11) – (D12) – (D13);
[black,thick] (D0) – (D1) – (D2) – (D3) – (D4) – (D5) – (D6) – (D7) – (D8) – (D9) – (D10) – (D11) – (D12) – (D13);

[black,very thick]
(S1_0) – (S1_1)
(S2_0) – (S2_1)
(S3_0) – (S3_1)
(S4_0) – (S4_1)

(S6_0) – (S6_1)
(S7_0) – (S7_1)
(S8_0) – (S8_1)
(S9_0) – (S9_1)

(S11_0) – (S11_1)
(S12_0) – (S12_1)
(S13_0) – (S13_1)
(S14_0) – (S14_1)
(S15_0) – (S15_1)
(S16_0) – (S16_1)
(S17_0) – (S17_1)
(S18_0) – (S18_1)
(S19_0) – (S19_1)

[class A] at (1.25,2.5) V;
[class A] at (3.125,3.75) V;
[class A] at (3.125,1.25) V;
[class A] at (1.25,6.25) V;
[class A] at (3.125,6.25) V;

[class B] at (7.5,7.5) S;

[class B] at (6.25,3.75) S;
[class B] at (6.25,3.75) S;
[class B] at (8.75,3.75) S;
[class B] at (8.75,1.25) S;
[class B] at (6.875,1.25) S;
[class B] at (1.25,9.375) S;
[class B] at (3.75,9.375) S;

After 17 Splits.

axis equal image,
scale only axis,
xlabel=$\vec{x}_1$, ylabel=$\vec{x}_2$,
xmin=0, xmax=10,
ymin=0, ymax=10,
xtick=, ytick=,

(D0)  at (-1,7.0);
(D1)  at (-0.1,7.45);
(D2)  at (1.6,8.2);
(D3)  at (4.0,8.5);
(D4)  at (4.6,8.1);
(D5)  at (4.75,7.3);
(D6)  at (4.65,6.3);
(D7)  at (4.4,4.7);
(D8)  at (4.5,3.8);
(D9)  at (4.7,2.6);
(D10) at (5.1,2.0);
(D11) at (5.5,1.3);
(D12) at (5.8,-0.2);
(D13) at (6.0,-1);

(S1_0) at (5,0);
(S1_1) at (5,10);

(S2_0) at (0,5);
(S2_1) at (5,5);

(S3_0) at (2.5,0);
(S3_1) at (2.5,5);

(S4_0) at (5,5);
(S4_1) at (10,5);

(S6_0) at (7.5,0);
(S6_1) at (7.5,5);

(S7_0) at (5,2.5);
(S7_1) at (10,2.5);

(S8_0) at (5,2.5);
(S8_1) at (10,2.5);

(S9_0) at (6.25,0);
(S9_1) at (6.25,2.5);

(S11_0) at (2.5,2.5);
(S11_1) at (5,2.5);

(S12_0) at (3.75,2.5);
(S12_1) at (3.75,5);

(S13_0) at (2.5,5);
(S13_1) at (2.5,10);

(S14_0) at (0,7.5);
(S14_1) at (2.5,7.5);

(S15_0) at (2.5,7.5);
(S15_1) at (5,7.5);

(S16_0) at (3.75,0);
(S16_1) at (3.75,2.5);

(S17_0) at (3.75,5);
(S17_1) at (3.75,7.5);

(S18_0) at (2.5,8.75);
(S18_1) at (5,8.75);

(S19_0) at (0,8.75);
(S19_1) at (2.5,8.75);

3 [surf,samples=30,domain=0:10,shader=interp] exp(-0.5*(((x-8)/5)^2 + ((y-8)/4)^2))/(sqrt(2*pi*5*4);

[gray,thin,step=1.25] (0,0) grid (10,10);

[side by side=Class2Class1] (D0) – (D1) – (D2) – (D3) – (D4) – (D5) – (D6) – (D7) – (D8) – (D9) – (D10) – (D11) – (D12) – (D13);
[black,thick] (D0) – (D1) – (D2) – (D3) – (D4) – (D5) – (D6) – (D7) – (D8) – (D9) – (D10) – (D11) – (D12) – (D13);

[black,very thick]
(S1_0) – (S1_1)
(S2_0) – (S2_1)
(S3_0) – (S3_1)
(S4_0) – (S4_1)

(S6_0) – (S6_1)
(S7_0) – (S7_1)
(S8_0) – (S8_1)
(S9_0) – (S9_1)

(S11_0) – (S11_1)
(S12_0) – (S12_1)
(S13_0) – (S13_1)
(S14_0) – (S14_1)
(S15_0) – (S15_1)
(S16_0) – (S16_1)
(S17_0) – (S17_1)
(S18_0) – (S18_1)
(S19_0) – (S19_1)

[class A] at (1.25,2.5)   V;
[class A] at (3.125,3.75) V;
[class A] at (3.125,1.25) V;
[class A] at (1.25,6.25)  V;
[class A] at (3.125,6.25) V;

[class B] at (7.5,7.5)   blackS;

[class B] at (6.25,3.75)  blackS;
[class B] at (6.25,3.75)  blackS;
[class B] at (8.75,3.75)  blackS;
[class B] at (8.75,1.25)  blackS;
[class B] at (6.875,1.25) blackS;
[class B] at (1.25,9.375) blackS;
[class B] at (3.75,9.375) blackS;

Probability Density Function.

[Computing Bounds on Probabilities]
Computing Bounds on Probabilities.
This figure illustrates the steps for computing bounds on \(p = \pterm\).
Our algorithm successively splits the input space to find regions that do not intersect the satisfaction boundary \(\parg\) ( line Class2Class1).
Green, orange, and grey rectangles (Class2/Class1/black!10!white) denote regions for which we could prove \(\parg\) (satisfaction) Class2, \(\fSAT*(\inp, \NN(\inp)) < 0\) (violation) Class1, or neither black!10!white, respectively.
By integrating the probability density \(\probdensity[\inp]\) in fig:probability-bounds-probability-mass (darker means higher density) over the green rectangles Class2, we obtain a lower bound on \(p\).
Similarly, we can integrate over the orange rectangles Class1 to construct an upper bound on \(p\).

Refining the input splitting from fig:probability-bounds-four-splits to fig:probability-bounds-nineteen-splits tightens the bounds on \(p\).

§.§ Probabilistic Verification Algorithm

<Ref> describes the  algorithm.

The centrepiece of  is the procedure  for computing bounds \(\ell_i^{(t)} \leq p_i \leq u_i^{(t)}\) on a probability \(p_i\) from <ref>.
Given \(\ell_i^{(t)} \leq p_i \leq u_i^{(t)}\), we apply interval arithmetic as introduced in <ref> to prove or disprove \(\fSAT(p_1, \ldots, p_v)\).
As described in <ref>, this analysis may be inconclusive.
In this case,  refines \(\ell_i^{(t)}, u_i^{(t)}\) to obtain \(\ell_i^{(t+1)}, u_i^{(t+1)}\) with \(\ell_i^{(t)} \leq \ell_i^{(t+1)} \leq p_i \leq u_i^{(t+1)} \leq u_i^{(t)}\).
We again apply interval arithmetic to \(\fSAT(p_1, \ldots, p_v)\), this time using \(\ell_i^{(t+1)}, u_i^{(t+1)}\).
If the result remains inconclusive, we iterate refining the bounds on each \(p_i\) until we obtain a conclusive result.
applies  for each \(p_i\) in parallel, making use of several CPU cores or several GPUs.
Our main contribution is the  algorithm for computing a converging sequence of lower and upper bounds on \(p_i\).

<Ref> describes  in detail.

Probabilistic Verification Problem as in <ref>, Batch Size \(N\)
(*[f]Launch \(v\) parallel instances of )$i \in \{1, \ldots, v\}$
${P\!B}_i \gets$ Launch $\ProbBounds(p_i, N)$

$t \in \Nats$
$i \in \{1, \ldots, v\}$
Gather $\ell_i^{(t)}, u_i^{(t)}$ from ${P\!B}_i$

\(\ell^{(t)}, u^{(t)} \gets \IntervalArithmetic(\fSAT, (\ell_1^{(t)}, u_1^{(t)}), \ldots, (\ell_v^{(t)}, u_v^{(t)}))\)
$\ell^{(t)} \geq 0$

$u^{(t)} < 0$

Probability $\pterm$, Batch Size $N$
$\mathrm{branches} \gets \{\mathcal{X}\}$*\(\mathcal{X}\) is the input space of \(\NN\)
$\ell^{(0)} \gets 0, u^{(0)} \gets 1$
$t \in \Nats$
$\mathrm{batch} \gets \Select(\mathrm{branches}, N)$ *e.g., branches \(\branch[i]\) with largest \(\prob[\inp](\branch[i])\)
$(\underline{\vec{y}}, \overline{\vec{y}}) \gets \ComputeBounds(\fSAT*(\cdot, \NN(\cdot)), \mathrm{batch})$
$(\mathrm{batch}, \Xsat[t], \Xviol[t]) \gets \Prune{}(\mathrm{batch}, \underline{\vec{y}}, \overline{\vec{y}})$
$\ell^{(t)} \gets \ell^{(t-1)} + \prob[\inp](\Xsat[t])$
$u^{(t)} \gets u^{(t-1)} - \prob[\inp](\Xviol[t])$
$\ell^{(t)}, u^{(t)}$ *Report new bounds to
$\mathrm{new} \gets \Split{}(\mathrm{batch})$ *e.g.,  <cit.>
$\mathrm{branches} \gets (\mathrm{branches} \setminus \mathrm{batch}) \cup \mathrm{new}$

§.§ Bounding Probabilities

Our  algorithm for deriving and refining bounds on a probability is described in detail in <ref> and illustrated in <ref>.
is a massively parallel input-splitting branch and bound procedure that leverages a bound propagation algorithm for non-probabilistic neural network verification ().
Since we only consider a single probability \(p_i\) as in <ref> in this section, we denote this probability as \(p = \pterm\).

receives \(p\) and a batch size \(N \in \Nats\) as input.
The algorithm iteratively computes \(\ell^{(t)}, u^{(t)} \in [0, 1]\), such that \(\ell^{(t)} \leq \ell^{(t')} \leq p \leq u^{(t')} \leq u^{(t)}\), \(\forall t, t' \in \Nats, t' \geq t\).

The following sections describe each step of  in detail.

§.§.§ Initialisation

Initially, we consider a single branch encompassing \(\NN\)'s entire input space \(\mathcal{X}\).
As in <ref>, we assume \(\mathcal{X}\) to be a bounded hyperrectangle.
We use the trivial bounds \(\ell^{(0)} = 0 \leq p \leq 1 = u^{(0)}\) as initial bounds on \(p\).

§.§.§ Selecting Branches

First, we select a batch of \(N \in \Nats\) branches.
In the spirit of XuZhangWangEtAl2021, we leverage the data parallelism of modern CPUs and GPUs to process several branches at once.
In iteration \(t=1\), the batch only contains the branch \(\mathcal{X}\).
Which branches we select determines how fast we obtain tight bounds on \(p\).
We propose two heuristics for selecting branches:

* : Inspired by  <cit.>, this heuristic selects the \(N\) branches \(\branch[i]\) with the largest \(\prob[\inp](\branch[i])\).
This heuristic is motivated by the observation that pruning these branches would lead to the largest improvement of \(\ell^{(t)}, u^{(t)}\).

* : This heuristic selects the \(N\) branches \(\branch[i]\) with the largest \(\frac{\prob[\inp](\branch[i])}{\log(e + \overline{\vec{y}}_i - \underline{\vec{y}}_i)}\) where \(\underline{\vec{y}}, \overline{\vec{y}}\) are as in .
The motivation for this heuristic is to select the branches with the largest probability while \(\underline{\vec{y}}, \overline{\vec{y}}\) are loose to focus on the most relevant regions of the input space, but greedily select branches that have tight bounds with the expectation that these will be pruned soon.

We compare  and  experimentally in <ref>.
The comparison reveals that  slightly speeds up  compared to .

§.§.§ Pruning

The next step is to prune those branches \(\branch[j] \in \mathrm{batch}\), for which we can determine that \(y = \parg\) is either certainly satisfied or certainly violated.
For this, we first compute \(\underline{\vec{y}} \leq \fSAT*(\cdot, \NN(\cdot)) \leq \overline{\vec{y}}\) for the entire \(\mathrm{batch}\) using a sound neural network verifier like interval arithmetic, *, or \(\alpha\)-CROWN <cit.>.

If \(\underline{\vec{y}}_j \geq 0\) (\(y \geq 0\) is certainly satisfied) or \(\overline{\vec{y}}_j < 0\) (\(y \geq 0\) is certainly violated), we can prune \(\branch[j]\) analogously to <ref>.
We collect the branches with \(\underline{\vec{y}}_j \geq 0\) in the set \(\Xsat[t]\) and the branches with \(\overline{\vec{y}}_j < 0\) in \(\Xviol[t]\), where \(t \in \Nats\) is the current iteration.

§.§.§ Updating Bounds

Let \(\Xsat*[t] = \bigcup_{t'=1}^{t} \Xsat[t']\) and \(\Xviol*[t] = \bigcup_{t'=1}^{t} \Xviol[t']\), where \(t \in \Nats\) is the current iteration.
Then, \(\ell^{(t)} = \prob[\inp](\Xsat*[t]) \leq p\).
Similarly, \(\ka^{(t)} = \prob[\inp](\Xviol*[t]) \leq \prob[\inp](\fSAT(\inp, \NN(\inp)) < 0) = 1 - p\).
Therefore, \(1 - \ka^{(t)} = u^{(t)} \geq p\).
Practically, we only have to maintain the current bounds \(\ell^{(t)}\) and \(u^{(t)}\) instead of the sets \(\Xsat*[t]\) and \(\Xviol*[t]\).

Because \(\Xsat*[t]\) and \(\Xviol*[t]\) are a union of disjoint hyperrectangles, exactly computing \(\prob[\inp](\Xsat*[t])\) and \(\prob[\inp](\Xviol*[t])\) is feasible for a large class of probability distributions, including discrete, uniform, and univariate continuous distributions, as well as Mixture Models and Bayesian Networks of such distributions, but not, for example, multivariate normal distributions.
Practically, one must also accommodate floating point errors in these computations.
While we do not address this in this paper, our approach can be extended to this end.

now reports the refined bounds \(\ell^{(t)} \leq p \leq u^{(t)}\) to , which applies interval arithmetic to determine whether \(\NN\) is feasible for the probabilistic verification problem.
If this remains inconclusive, we proceed with splitting the selected branches.

§.§.§ Splitting

Splitting refines a branch \(\branch = \HR{\inp}\) by selecting a dimension \(d \in \{1, \ldots, n\}\) to split.
We first describe how \(d\) is split before discussing how to select \(d\).
A dimension can encode several types of variables.
We consider continuous variables, such as normalised pixel values, integer variables, such as age, and dimensions containing one indicator of a one-hot encoded categorical variable like gender.
The type of variable encoded in \(d\) determines how we split \(d\).

* For a continuous variable, we bisect \(\branch\) along \(d\) resulting in two new branches \(\HR[']{\inp}\) and \(\HR['']{\inp}\).
Concretely, \(\underline{\inp}'_{d'} = \underline{\inp}''_{d'} = \underline{\inp}_{d'}\) and \(\overline{\inp}'_{d'} = \overline{\inp}''_{d'} = \overline{\inp}_{d'}\) for all \(d' \in \{1, \ldots, n\} \setminus \{d\}\) while \(\overline{\inp}'_d = \underline{\inp}''_d = (\underline{\inp}_d + \overline{\inp}_d)/2\), \(\underline{\inp}'_d = \underline{\inp}_d\), and \(\overline{\inp}''_d = \overline{\inp}_d\).

* For integer variables, we bisect along \(d\) to obtain \(\HR[']{\inp}\), \(\HR['']{\inp}\) and round \(\overline{\inp}'_d\) to the next smaller integer while rounding \(\underline{\inp}''_d\) to the next larger integer.

* For a one-hot encoded categorical variable \(A\) encoded in the dimensions \(D \in \{1, \ldots, n\}\) with \(d \in D\), we create one split where \(A\) is equal to the category represented by \(d\) and one where \(A\) is different from this category.
Formally, \(\underline{\vec{x}}_d' = \overline{\vec{x}}_d' = 1\) and \(\underline{\vec{x}}_{d'}' = \overline{\vec{x}}_{d'}' = 0\) for \(d' \in D \setminus \{d\}\) defines \(\HR[']{\vec{x}}\).

For \(\HR['']{\vec{x}}\), we set \(\underline{\vec{x}}_{d}'' = \overline{\vec{x}}_{d}'' = 0\) and leave the remaining values are they are in \(\underline{\vec{x}}\) and \(\overline{\vec{x}}\)[This splitting procedure eventually creates a new branch where all dimensions are set to zero. This branch has zero probability and can be discarded immediately.].

In any case, we need to ensure not to select \(d\) if \(\underline{\inp}_d = \overline{\inp}_d\).
We introduce two heuristics for selecting a dimension:

* : This well-known heuristic <cit.> selects the dimension with the largest edge length \(\overline{\inp}_d - \underline{\inp}_d\).

* : We use a variant of the  heuristic of BunelLuTurkaslanEtAl2020.
The idea of  is to estimate the improvement in bounds that splitting dimension \(d\) yields by using a yet less expensive technique than .
Our variant of  uses interval arithmetic, assuming that we use  or \(\alpha\)-CROWN for .
Let \(\HR[^{(d, 1)}]{\inp}\) and \(\HR[^{(d, 2)}]{\inp}\) be the two new branches originating from splitting dimension \(d \in \{1, \ldots, n\}\) and let \(\underline{y}^{(d, 1)}, \overline{y}^{(d, 1)}\), \(\underline{y}^{(d, 2)}, \overline{y}^{(d, 2)}\) be the bounds that interval arithmetic computes on \(\fSAT*(\cdot, \NN(\cdot))\) for these branches.
Our  selects \(d = \argmax_{d \in \{1, \ldots, n\}} \tilde{y}^{(d)}\), where \(\tilde{y}^{(d)} = \max(\max(\underline{y}^{(d, 1)}, \underline{y}^{(d, 2)}), -{\min(\overline{y}^{(d, 1)}, \overline{y}^{(d, 2)})})\).

While  is more theoretically accessible,  is practically advantageous, as discussed in <ref>.

§ THEORETICAL ANALYSIS

In this section, we prove that  is a sound probabilistic verification algorithm when instantiated with a suitable  procedure.
We also prove that  instantiated with , , and interval arithmetic for  is complete under mild assumptions on the probabilistic verification problem.
Soundness and completeness are defined in <ref>.
As in <ref>, we omit index and superscripts when considering only a single probability \(p = \pterm\).

§.§ Soundness

First, we prove that  produces sound bounds on \(p\) when using a sound  procedure, that is, a procedure that computes valid bounds, such as interval arithmetic or *.
The soundness of  then follows immediately.

Let \(N \in \Nats\) be a batch size and assume  produces valid bounds.
Let \({\left((\ell^{(t)}, u^{(t)})\right)}_{t \in \Nats}\) be the iterates of \(\ProbBounds(\pterm, N)\).
It holds that \(\ell^{(t)} \leq \pterm \leq u^{(t)}\) for all $t\in\mathbb{N}$.

Let \(t \in \Nats\) and let \(\Xsat[t]\) and \(\Xviol[t]\) be as in <ref>.
computes \(\ell^{(t)}\) as the total probability of all previously pruned satisfied branches \(\Xsat*[t] = \bigcup_{t'=1}^t \Xsat[t']\).
Similarly, \(u^{(t)} = 1 - \ka^{(t)}\) where \(\ka^{(t)}\) is the total probability of all previously pruned violated branches \(\Xviol*[t] = \bigcup_{t'=1}^t \Xviol[t']\).
Since we assumed that  produces valid bounds,  only prunes branches that are actually satisfied or violated.

Therefore, \(\Xsat*[t] \subseteq \{\inp \in \mathcal{X} \mid \parg\}\) and \(\Xviol*[t] \subseteq \{\inp \in \mathcal{X} \mid \fSAT*(\inp, \NN(\inp)) < 0\}\).
From this, it follows directly that

\begin{align*}
\ell^{(t)} &= \prob[\inp](\Xsat*[t])
\leq \pterm \\
\ka^{(t)} &= \prob[\inp](\Xviol*[t])
\leq \prob[\inp](\fSAT*(\inp, \NN(\inp)) < 0),
\end{align*}

which implies \(u^{(t)} = 1 - \ka^{(t)} \geq 1 - \prob[\inp](\fSAT*(\inp, \NN(\inp)) < 0) = \pterm\).
This shows that  is sound.

is sound when using a  procedure that computes valid bounds.

<Ref> follows from <ref> and the soundness of interval arithmetic <cit.>[We include relevant theorems from MooreKearfottCloud2009 in <ref> for reference.].

§.§ Completeness

We study the completeness of .
Concretely, we prove that  instantiated with , interval arithmetic for , and  is complete under a mildly restrictive condition on <ref>.

Let \(v\), \(\fSAT\), \(\fSAT*[i]\), and \(\inp[i]\) be as in <ref>.
Assume \(\fSAT(p_1, \ldots, p_v) \neq 0\) and \(\forall i \in \{1, \ldots, v\}: \prob[\inp[i]](\pfun[i]{} = 0) = 0\).

To prove the completeness of , we first establish that  produces a sequence of lower and upper bounds that converge towards each other.
Intuitively, we require <ref> since converging bounds on \(\fSAT(p_1, \ldots, p_n)\) are insufficient for proving \(\fSAT(p_1, \ldots, p_n) \geq 0\) if \(\fSAT(p_1, \ldots, p_n) = 0\) <cit.>.
However, excluding \(\fSAT(p_1, \ldots, p_v) = 0\) is only mildly restrictive as we can always tighten the constraint to \(\fSAT(p_1, \ldots, p_n) \geq \varepsilon\) for an arbitrarily small \(\varepsilon > 0\).
The second assumption \(\prob[\inp](\pfun{} = 0) = 0\) is only mildly restrictive for similar reasons.
In particular, we can tighten \(\pterm\) to \(\prob[\inp](\pfun{} \geq \varepsilon)\) for \(\varepsilon > 0\) such that \(\prob[\inp](\pfun{} = \varepsilon) = 0\).
Such a \(\varepsilon > 0\) exists because \(\prob[\inp](\pfun{} = 0) = 0\) means that the satisfaction boundary has positive volume, but any neural network has only finitely many flat regions that can produce a satisfaction boundary of positive volume.

Let \(N \in \Nats\) be a batch size.
Let \({\left((\ell^{(t)}, u^{(t)})\right)}_{t \in \Nats}\) be the iterates of \(\ProbBounds{}(\pterm, N)\) instantiated with , interval arithmetic for  and .
Assume \(\prob[\inp](\fSAT*(\inp, \NN(\inp)) = 0) = 0\) as in <ref>.

\begin{equation*}
\lim_{t \to \infty} \ell^{(t)} = \lim_{t \to \infty} u^{(t)} = \pterm.
\end{equation*}

Purely Discrete Input Spaces.
In the following, we assume that the input space \(\mathcal{X}\) of \(\NN\) contains at least one continuous variable.
Otherwise, \(\mathcal{X}\) contains only finitely many discrete values.
This implies that  eventually reaches an iteration where no branch can be split further.
However, in this iteration, all branches are points in the input space, for which interval arithmetic computes \(\underline{\vec{y}} = \overline{\vec{y}}\) <cit.>.
In turn, this implies that all branches are pruned, making  complete.

We require the following intermediate result for proving <ref>.
We write \(\branch \leadsto \branch[t]\) if \(\branch[t]\) is a branch in iteration \(t\) of  that originates from splitting \(\branch\), meaning that \(\branch[t] \subset \branch\).

Let \(N \in \Nats\).
instantiated as in <ref> satisfies

\begin{equation*}
\forall t \in \Nats:
\forall \branch \in \mathrm{branches}^{(t)}:
\prob[\inp](\branch) > 0 \Longrightarrow
\exists t' \geq t: \branch \in \Prob{}(\mathrm{branches}^{(t')}, N),
\end{equation*}

where \(\mathrm{branches}^{(t)}\) is the value of the \(\mathrm{branches}\) variable of  in iteration \(t\).

Let \(N \in \Nats\) and let \(\branch \notin \Prob{}(\mathrm{branches}^{(t)}, N)\) in iteration \(t \in \Nats\) of  with \(\prob[\inp](\branch) > 0\).

There are at least \(N\) branches \(\branch[t]'\) in iteration \(t\) with \(\prob[\inp](\branch[t]') \geq \prob[\inp](\branch)\).
We show

\begin{equation}
\exists t' > t:
\forall \branch[t]', \prob[\inp](\branch[t]') \geq \prob[\inp](\branch):
\underbrace{%
\forall \branch[t']', \branch[t]' \leadsto \branch[t']':
\prob[\inp](\branch[t']') < \prob[\inp](\branch)
\label{eqn:evtl-prob-low}
\end{equation}

Let \(\branch[t]'\) be a branch in iteration \(t\) with \(\prob[\inp](\branch[t]') \geq \prob[\inp](\branch[t])\).
We first show that there is an iteration \(t' > t\) such that \((\ast)\) holds for \(\branch[t]'\).

First of all, if \(\branch[t]'\) is pruned by  in iteration \(t\), then there are no new branches originating from \(\branch[t]'\), so that \((\ast)\) holds vacuously.
Otherwise,  splits \(\branch[t]'\).

Without loss of generality, assume that the dimension selected for splitting encodes a continuous variable.
This does not harm generality since discrete variables in a bounded input space can only be split finitely often and will, therefore, eventually become unavailable for splitting.
Since we split continuous variables by bisection, we have that the volume of all branches \(\branch[t']'\) originating from \(\branch[t]'\) decreases towards zero as \(t'\) increases.

As stated in <ref>, we assume that all continuous random variables admit a probability density function.
This implies that the probability in all branches \(\branch[t']'\) originating from \(\branch[t]'\) decreases towards zero as the volume decreases towards zero.
Therefore, there is a \(t'\), such that \((\ast)\) is satisfied for \(\branch[t]'\).
Since there are only finitely many branches in any iteration of , the above implies that <ref> is satisfied.
In turn, this directly implies that \(\branch\) is eventually selected by , proving <ref>.

We first prove that \(\lim_{t \to \infty} \ell^{(t)} = \pterm\).
The convergence of the upper bound, \(\lim_{t \to \infty} u^{(t)} = \pterm\) follows with a similar argument.

Let \(\Xsat<\ast> = \{\inp \in \mathcal{X} \mid \pfun > 0\}\), where \(\mathcal{X}\) is the input space of \(\NN\).
Note that \(\prob[\inp](\Xsat<\ast>) = \prob[\inp](\{\inp \in \mathcal{X} \mid \parg\})\) due to <ref>.
Further, let \(\Xsat*[t]\) be as in the proof of <ref> and recall \(\ell^{(t)} = \prob[\inp](\Xsat*[t])\).

First, we give an argument why the limit \(\lim_{t \to \infty} \ell^{(t)}\) exists.
Due to <ref>, \(\ell^{(t)}\) is bounded from above.
Furthermore, \(\ell^{(t)}\) is non-decreasing in \(t\) since  only adds elements to \(\Xsat*[t]\).
Therefore, \(\lim_{t \to \infty} \ell^{(t)}\) exists.

Given this, we now show

[](*[t])        [](<∗>)

[](*[t]) - [](<∗>)        0

[](<∗> ∖*[t])        0.

We prove <ref> by showing that \(\Xsat**[t] = \Xsat<\ast> \setminus \Xsat*[t]\) has diminishing volume.
Since we assumed the input space to contain at least one continuous variable and assumed all continuous random variables to admit a density function, diminishing volume implies diminishing probability.

With the goal of obtaining a contradiction, assume \(\lim_{t \to \infty} \vol(\Xsat**[t]) > 0\), where \(\vol\) denotes the volume of a set.

Let \(t \in \Nats\).
Since the branches maintained by  form a partition of the input space, there is a branch \(\branch[t]\) in iteration \(t\) of  such that \(\Xsat** \cap \branch[t] \neq \emptyset\).

A central argument in the following is that  eventually splits every dimension since splitting decreases the edge length \(\overline{\vec{x}}_d - \underline{\vec{x}}_d\) of the dimension \(d\) being split.
Due to this and due to <ref>, we eventually obtain \(\branch[t']\) in iteration \(t' > t\) with \(\branch[t'] \subseteq \Xsat** \cap \branch[t]\).
Also, since  eventually splits every dimension, Theorem 6.1 of MooreKearfottCloud2009 applies, which states that the bounds \(\underline{y} \leq \pfun \leq \overline{y}\) produced by interval arithmetic converge towards \(\pfun\).
However, this implies that there is an iteration \(t'' > t'\) in which  considers a branch \(\branch[t''] \subset \Xsat**\) for which \(\underline{y} > 0\), which means that  prunes \(\branch[t'']\).
This contradicts the construction of\(\Xsat**\).

With this contradiction, we have shown \(\lim_{t \to \infty} \ell^{(t)} = \pterm\).
Convergence of the upper bound \(u^{(t)}\) follows from an analogous argument on \(\ka^{(t)}\) as in the proof of <ref> where \(u^{(t)} = 1 - \ka^{(t)}\).
This establishes <ref>.

When instantiated with  as in <ref>,  is complete for verification problems satisfying <ref>.

Let \(\NN\), \(\fSAT\), \(\fSAT*[1], \ldots, \fSAT*[v]\) be as in <ref>.

First, consider \(\fSAT(p_1, \ldots, p_v) > 0\).
As a consequence of Theorem 6.1 of MooreKearfottCloud2009 and <ref>, the bounds \(\ell \leq \fSAT(p_1, \ldots, p_v) \leq u\) produced by interval arithmetic converge towards \(\fSAT(p_1, \ldots, p_v)\).
This implies that eventually \(\ell > 0\), meaning that  eventually proves \(\fSAT(p_1, \ldots, p_v) \geq 0\).

If \(\fSAT(p_1, \ldots, p_v) < 0\) we eventually obtain \(u < 0\) with the same argument as above.
Since \(u < 0\) disproves \(\fSAT(p_1, \ldots, p_v) \geq 0\),  is complete for probabilistic verification problems satisfying <ref>.

§ EXPERIMENTS

We apply  (Algorithm <ref>) to verify the demographic parity fairness notion, count the number of safety violations of a neural network controller in a safety-critical system, and quantify the robustness of a neural network.
We use the  and  heuristics and use * for .
<Ref> gives an overview of the benchmarks considered in this section.
All verification problems are defined formally in <ref>.

As our results show,  outpaces the probabilistic verification algorithms  <cit.> and  <cit.>.
Additionally, we show that  compares favourably to the  <cit.> and  <cit.> algorithms for #DNN verification <cit.>, which corresponds to probabilistic verification with uniformly distributed inputs.

While no code is publicly available for , running  is very computationally expensive.
To enable a faithful comparison, we run our experiments on less powerful hardware (HW1) compared to the hardware used by ConverseFilieriGopinathEtAl2020,MarzariRoncolatoFarinelli2023 and compare the runtime of our algorithms to the runtimes reported by these authors.
To obtain comparable runtimes across the experiments in this section, we run all experiments on HW1.
The results of  and  were similarly obtained on HW1.
Additional experiments on more powerful hardware are contained in <ref>

To test the limits of , we introduce a new, challenging benchmark:
MiniACSIncome is based on the ACSIncome dataset <cit.> and consists of datasets of varying input dimensionality, probability distributions for these datasets, and neural networks of varying sizes.
solves seven of eight instances in MiniACSIncome within an hour.

Benchmark     Input Dimension     Input Distributions     Network Size     Source
2*FairSquare        2*2–3     independent             2*1$\times$1, 1$\times$2     2*<cit.>

2 Bayesian Networks

ACAS Xu           5        uniform                 6$\times$50     <cit.>

2*MiniACSIncome     2*1–8      2*Bayesian Network     1$\times$10–10000     2*Own

1$\times$10 – 10$\times$10

Benchmarks. Network size is the size of the neural network given as #layers$\times$layer size.

Hardware and Implementation.

We implement  in Python, leveraging PyTorch <cit.> and auto_LiRPA <cit.>.
Our code and benchmarks are available at <https://github.com/sen-uni-kn/probspecs>.
We run all experiments on a desktop running Ubuntu 22.04 with an Intel i7–4820K CPU, 32 GB of memory, and no GPU (HW1).
<Ref> compares our hardware to the hardware used by MarzariRoncolatoFarinelli2023,ConverseFilieriGopinathEtAl2020.

§.§ FairSquare Benchmark

AlbarghouthiDAntoniDrewsEtAl2017 evaluate their  algorithm on an application derived from the Adult <cit.> dataset.
In particular, they verify the fairness of three small neural networks with respect to a person's sex under three different distributions of the network input: a distribution of entirely independent univariate variables and two Bayesian Networks.
We replicate these experiments to compare  to .

The FairSquare benchmark assumes that all input variables except sex are unbounded, although some variables, like a person's age, have natural bounds.
To apply , we set the bounds of the input space such that the cumulative distribution functions of the input variables are $0$ at the lower bound and $1$ at the upper bound, up to machine precision.

<Ref> contains the results of our experiment.
Our verification results match the results of .
Regarding runtime,  significantly outperforms .
In particular,  verifies NN3,2 in several seconds, while  requires more than seven minutes for the independent input distribution and exceeds the time budget of 15 minutes for the Bayesian network input distributions.

A comparison with  on an extended set of benchmarks is contained in <ref>.

2cBenchmark Instance
2cRuntime (s)

$\NN$     $\prob[\inp]$
(Ours)      <cit.>     Fair?
NN2,1     independent     1.3     2.1

NN2,1     Bayes Net 1     2.2     75.5

NN2,1     Bayes Net 2     3.6     50.5

NN2,2     independent     1.8     4.2

NN2,2     Bayes Net 1     3.3     33.3

NN2,2     Bayes Net 2     4.9     58.6

NN3,2     independent     2.2     450.5

NN3,2     Bayes Net 1     4.4     TO

NN3,2     Bayes Net 2     6.9     TO

FairSquare Benchmark.
NN\(i\),\(j\) denotes a neural network with \(i\) inputs and a single hidden layer of \(j\) neurons.
TO indicates timeout (900s).

§.§ ACAS Xu Safety

The ACAS Xu networks <cit.> are a suite of 45 networks, together forming a collision avoidance system for crewless aircraft.
We reproduce the ACAS Xu #DNN verification <cit.> experiments of MarzariRoncolatoFarinelli2023.
In particular, we seek to quantify the number of violations (violation rate) of the ACAS Xu networks that violate the property \(\phi_2\) introduced by KatzBarrettDillEtAl2017.
This corresponds to computing bounds on <ref> under a uniform distribution of \(\inp\).

We compare  (<ref>) to the  and  algorithms for #DNN verification.
<cit.> computes the violation rate exactly, while  <cit.> computes an upper bound on the violation rate that is sound with a certain probability.
In contrast,  provides sound bounds on the violation rate at any time during its execution.
Therefore, we can run  only for a particular time budget or abort  when the bounds become sufficiently tight.

<Ref> compares  (Algorithm <ref>) with  and  for the ACAS Xu property \(\phi_2\) <cit.> and the networks investigated by MarzariRoncolatoFarinelli2023.

For all three networks,  can tighten the bounds to a margin of less than $0.007$ within one hour, while  requires at least four hours to compute the exact violation rate.
In comparison to ,  produces tighter sound upper bounds within 10 seconds in two of three cases, while  requires at least 57 seconds to derive a probably sound upper bound for these cases.
The extended comparison in <ref> reveals that, in 13 from a total of 36 cases,  computes tighter sound bounds faster than  computes a probably sound upper bound.

3c (Ours)
2c <cit.>
2c <cit.>

10s     1m     1h
2c99.9% confid.

\(\NN\)
\(\ell, u\)     \(\ell, u\)     \(\ell, u\)
VR     Rt     \(u\)     Rt

$N_{4,3}$     $  0.19\%,   2.90\%$     $  0.62\%,   2.28\%$     $  1.12\%,   1.75\%$     $1.43\%$     8h 46m      $3.61\%$     64.5s

$N_{4,9}$     $  0.00\%,   3.39\%$     $  0.00\%,   1.65\%$     $  0.08\%,   0.30\%$     $0.15\%$     12h 21m     $0.73\%$     19.5s

$N_{5,8}$     $  0.95\%,   4.12\%$     $  1.58\%,   3.11\%$     $  1.97\%,   2.57\%$     $2.20\%$      4h 35m      $4.52\%$     57.0s

Comparison of , , and .
We run  with different time budgets (10s, 1m, 1h) and report the lower and upper bounds (\(\ell, u\)) computed within this time budget.
In contrast,  computes the exact probabilities (VR), and  computes a 99.9% confidence (confid.) upper bound.
The runtimes (Rt) of  are taken from MarzariRoncolatoFarinelli2023.

§.§ ACAS Xu Robustness

We replicate the experiments of ConverseFilieriGopinathEtAl2020 who apply  to quantify the robustness of ACAS Xu network \(N_{1,1}\) <cit.> under adversarial perturbations.
As ConverseFilieriGopinathEtAl2020, we consider 25 reference inputs — five for each class — and allow these reference inputs to be perturbed in the first two dimensions by at most \(5\%\) of the diameter of the input space in the respective dimension.
To compute bounds on the output distribution, we bound the probability of each of the five ACAS Xu classes for each of the 25 inputs.
Since the ACAS Xu training data is not publicly available, we sample the reference inputs randomly.

The mean runtime of  for these 125 verification problems is 18 seconds (median: 5.1s, maximum: 173.4s).

Therefore,  substantially outperforms , for which ConverseFilieriGopinathEtAl2020 report an average runtime of 33 minutes per verification problem while running their experiments on superior hardware.
<Ref> contains the full results.

§.§ MiniACSIncome

To test the limits of , we introduce the MiniACSIncome benchmark.
MiniACSIncome is derived from the ACSIncome <cit.> dataset, a replacement of the Adult dataset <cit.> that is better suited for fair machine learning research.
ACSIncome is based on US census data.
The task is to predict whether a person's yearly income exceeds $50 000 using features such as the person's age, sex, weekly work hours, education, etc.
Our benchmark provides probabilistic verification problems of various degrees of difficulty.
We apply  to MiniACSIncome and show that it outperforms a naïve approach for solving MiniACSIncome.

To create probabilistic verification problems of increasing difficulty, we consider an increasing number of input variables from ACSIncome.
The smallest instance, MiniACSIncome-1, only contains the binary SEX variable. In contrast, the largest instance, MiniACSIncome-8, contains SEX and seven more variables from ACSIncome, including age, education, and working hours per week.
We train a neural network with a single layer of ten neurons for each MiniACSIncome-$i$, \(i \in \{1, \ldots, 8\}\).

For MiniACSIncome-4, we additionally train deeper and wider networks to investigate the scalability of  with respect to the network size.

We fit a Bayesian Network to the MiniACSIncome-8 dataset to obtain an input distribution.
We use this distribution for all MiniACSIncome-$i$ instances by taking the variables not contained in MiniACSIncome-$i$ as latent variables.
The verification problem is then to verify the demographic parity of a neural network with respect to SEX under this input distribution.
<Ref> contains additional details on the MiniACSIncome benchmark.

Naïve Approach.

Since all variables in ACSIncome are discrete, we can theoretically verify the demographic parity of a neural network by enumerating all discrete values in the input space of MiniACSIncome-$i$.
However, as expected, the number of discrete values scales exponentially with the number of input variables.
While MiniACSIncome-1 contains only two values, MiniACSIncome-4 already contains 38 000  values, and MiniACSIncome-8 contains 2.4 billion values.
As we demonstrate next,  can explore this space more efficiently than the naïve enumeration approach.

<Ref> (left) displays the runtime of  and the naïve enumeration approach for MiniACSIncome-1 to MiniACSIncome-8.
While enumeration is faster than  for input spaces of up to three variables, where the network can be evaluated for all discrete values in one batch, enumeration falls behind  as soon as this becomes infeasible.
can solve MiniACSIncome for up to seven input variables in less than 30 minutes, only exceeding the timeout of one hour for MiniACSIncome-8.

Effect of Network Size.

<Ref> (right) displays the runtime of  for MiniACSIncome-4 networks of various sizes.
The studied networks include wide single-layer networks of up to 10 000 neurons and deep networks of up to 10 layers of 10 neurons.
As the figure shows,  is largely unaffected by the size of the MiniACSIncome-4 networks.
This is unexpected since the network size indirectly determines the performance of  through the complexity of the decision boundary.
However, larger networks need not necessarily have a more complex decision boundary, and large networks do not provide a performance benefit for MiniACSIncome-4, as discussed in <ref>.
Thoroughly exploring the impacts of network size requires more intricate datasets for which larger networks actually provide a benefit.

title=Scalability in the Input Size,
xmin=1, xmax=8,
ymin=0, ymax=3600,
xlabel=# Input Variables, ylabel=Runtime (min),
ylabel shift=-25pt,
legend pos=north west,
legend cell align=left,
[Colors-B,very thick,mark=x] table[x=InputVariables,y=EnumerateRuntime,col sep=comma] data/MiniACSIncomeInputSize.csv;

[Colors-A,very thick,mark=x] table[x=InputVariables,y=VerifyRuntime,col sep=comma] data/MiniACSIncomeInputSize.csv;

title=Scalability in the Network Size,

ymin=0, ymax=3600,
scaled ticks=false,

xlabel=# Parameters, ylabel=,
ylabel shift=-25pt,
legend pos=north west,
legend cell align=left,
[Colors-A,very thick,mark=x] table[x=NumParameters,y=Runtime,col sep=comma] data/MiniACSIncomeNumParameters.csv;

MiniACSIncome Results.
The left plot depicts the runtime of  and the naïve enumeration approach on MiniACSIncome-1 to MiniACSIncome-8 for a network of \(10\) neurons.
The right plot depicts the runtime of  for MiniACSIncome-4 networks of varying sizes.

§ CONCLUSION AND FUTURE WORK

As <ref> shows, our  algorithm for the probabilistic verification of neural networks significantly outpaces existing algorithms for probabilistic verification.
We obtain this speedup by recombining existing techniques for non-probabilistic neural network verification with the ideas from the  <cit.> fairness verifier.
Our MiniACSIncome benchmark provides a testbed for future probabilistic verification algorithms.

A promising direction for such verifiers is neuron branching <cit.>, which fuels the current state-of-the-art non-probabilistic neural network verifiers.
However, using neuron branching for probabilistic verification is challenging because it creates non-hyperrectangular regions in the input space that are expensive to integrate over.
Alternative directions include new branching and splitting heuristics, as well as switching to successively more precise bound propagation techniques during the course of .

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

howpublishedUCI Machine Learning Repository

titleOptuna: A Next-generation Hyperparameter Optimization Framework

journaltitleProc. ACM Program. Lang.
titleFairSquare: probabilistic verification of program fairness

journaltitleMath. Program.
titleStrong mixed-integer programming formulations for trained neural networks

booktitleCAV (1)
seriesLecture Notes in Computer Science
titleImproved Geometric Path Enumeration for Verifying ReLU Neural Networks

titleQuantitative Verification of Neural Networks and Its Security Applications

MIT Press

titleFairness and Machine Learning

journaltitleProc. ACM Program. Lang.
titleProbabilistic Verification of Fairness Properties via Concentration

seriesLecture Notes in Computer Science
titleSpecRepair: Counter-Example Guided Safety Repair of Deep Neural Networks

AAAI Press

titleProbabilistic Inference in Hybrid Domains by Weighted Model Integration

Computer Vision Foundation / IEEE Computer Society

titleAnalytic Expressions for Probabilistic Moments of PL-DNN With Gaussian Input

titleFairify: Fairness Verification of Neural Networks

seriesCEUR Workshop Proceedings
titleProvable Fairness for Neural Network Models Using Formal Verification

titleThe Fourth International Verification of Neural Networks Competition (VNN-COMP 2023): Summary and Results

journaltitleJ. Mach. Learn. Res.
titleBranch and Bound for Piecewise Linear Neural Network Verification

AUAI Press

seriesProceedings of Machine Learning Research
titleLagrangian Decomposition for Neural Network Verification

IEEE Computer Society

booktitleIEEE Symposium on Security and Privacy Workshops
titleAudio Adversarial Examples: Targeted Attacks on Speech-to-Text

seriesLecture Notes in Computer Science
titleMaximum Resilience of Artificial Neural Networks

titleProbabilistic Symbolic Analysis of Neural Networks

titleRetiring Adult: New Datasets for Fair Machine Learning

titleA DPLL(T) Framework for Verifying Deep Neural Networks

titleFairness through awareness

Association for Computational Linguistics

booktitleACL (2)
titleHotFlip: White-Box Adversarial Examples for Text Classification
family=European Parliament,

titleEU AI Act: first regulation on artificial intelligence

titleCertifying and Removing Disparate Impact

titleComplete Verification via Multi-Neuron Relaxation Guided Branch-and-Bound

booktitleESEC/SIGSOFT FSE
titleFairness testing: testing software for discrimination

IEEE Computer Society

booktitleIEEE Symposium on Security and Privacy
titleAI2: Safety and Robustness Certification of Neural Networks with Abstract Interpretation

titleBenchmarking Neural Network Robustness to Common Corruptions and Perturbations

Computer Vision Foundation / IEEE

titleNatural Adversarial Examples

titleDEEPSPLIT: An Efficient Splitting Method for Neural Network Verification via Indirect Effect Analysis

titleGoogle's Cloud Vision API is Not Robust to Noise

titleTowards Proving the Adversarial Robustness of Deep Neural Networks

booktitleCAV (1)
seriesLecture Notes in Computer Science
titleReluplex: An Efficient SMT Solver for Verifying Deep Neural Networks

titleAdam: A Method for Stochastic Optimization

titleCounterfactual Fairness

booktitle50 Years of Integer Programming 1958-2008 - From the Early Years to the State-of-the-Art
noteFirst published in Econometrica 28 (1960)
titleAn Automatic Method for Solving Discrete Programming Problems

titleThe #DNN-Verification Problem: Counting Unsafe Inputs for Deep Neural Networks

AAAI Press

titleEnumerating Safe Regions in Deep Neural Networks with Provable Probabilistic Guarantees

titleScaling #DNN-Verification Tools with Efficient Bound Propagation and Parallel Computing

titleFETA: Fairness Enforced Verifying, Training, and Predicting Algorithms for Neural Networks

titleIntroduction to Interval Analysis

titleA Unified Framework for Probabilistic Verification of AI Systems via Weighted Model Integration

journaltitleDiscret. Optim.
titleBranch-and-bound algorithms: A survey of recent advances in searching, branching, and pruning

titleThe Third International Verification of Neural Networks Competition (VNN-COMP 2022): Summary and Results

journaltitleProc. ACM Program. Lang.
titlePRIMA: general and precise neural network certification via scalable convex hull approximations

titleScaling the Convex Barrier with Active Sets

titlePyTorch: An Imperative Style, High-Performance Deep Learning Library

titleDiscrimination-aware data mining

seriesLecture Notes in Computer Science
titleAn Abstraction-Refinement Approach to Verification of Artificial Neural Networks

titleLearning Certified Individually Fair Representations

journaltitleProc. ACM Program. Lang.
titleAn Abstract Domain for Certifying Neural Networks

booktitleICLR (Poster)
titleIntriguing properties of neural networks

booktitleICLR (Poster)
titleEvaluating Robustness of Neural Networks with Mixed Integer Programming

booktitleCAV (1)
seriesLecture Notes in Computer Science
titleVerification of Deep Convolutional Neural Networks Using ImageStars

booktitleCAV (1)
seriesLecture Notes in Computer Science
titleNNV: The Neural Network Verification Tool for Deep Neural Networks and Learning-Enabled Cyber-Physical Systems

journaltitleProc. ACM Program. Lang.
titlePerfectly parallel fairness certification of neural networks

titleEfficient Formal Safety Analysis of Neural Networks

USENIX Association

booktitleUSENIX Security Symposium
titleFormal Security Analysis of Neural Networks using Symbolic Intervals

seriesProceedings of Machine Learning Research
titlePROVEN: Verifying Robustness of Neural Networks with a Probabilistic Approach

seriesProceedings of Machine Learning Research
titleTowards Fast Computation of Certified Robustness for ReLU Networks

seriesProceedings of Machine Learning Research
titleProvable Defenses against Adversarial Examples via the Convex Outer Adversarial Polytope

titleMarabou 2.0: A Versatile Formal Analyzer of Neural Networks

titleAutomatic Perturbation Analysis for Scalable Certified Robustness and Beyond

titleFast and Complete: Enabling Complete Neural Network Verification with Rapid and Massively Parallel Incomplete Verifiers

titleGeneral Cutting Planes for Bound-Propagation-Based Neural Network Verification

titleEfficient Neural Network Robustness Certification with General Activation Functions

§ PROBABILISTIC VERIFICATION PROBLEMS

This sections contains the formal definitions of all probabilistic verification problems in this paper.
First of all, we provide a detailed derivation of <ref>

[Derivation of <ref>]
Let \(\mathcal{X} \subset \Reals^n\), \(A \subset \{1, \ldots, n\}\), \(a \in A\), \(\NN: \Reals^n \to \Reals^2\) and \(\gamma\) be as in <ref>.

We rewrite <ref> as

[x⃗](() =  | is disadvantaged)
[x⃗](() =  | is advantaged)
≥    γ

[](()_1 - ()_2 ≥0 |_a ≤0)
[](()_1 - ()_2 ≥0 |_a ≥1)
≥    γ

[](()_1 - ()_2 ≥0 ∧_a ≤0) / [](_a ≤0)
[](()_1 - ()_2 ≥0 ∧_a ≥1) / [](_a ≥1)
≥    γ

[](min(()_1 - ()_2, -_a) ≥0) / [](-_a ≥0)
[](min(()_1 - ()_2, _a - 1) ≥0) / [](_a - 1≥0)
≥    γ

[](*[1](, ()) ≥0) / [](*[2](, ()) ≥0)
[](*[3](, ()) ≥0) / [](*[4](, ()) ≥0)
- γ   ≥    0

⟺    (
[](*[1](, ()) ≥0),
[](*[4](, ()) ≥0)
)    ≥    0

where \(\fSAT\), \(\fSAT*[1]\), \(\fSAT*[2]\), \(\fSAT*[3]\), and \(\fSAT*[4]\) are as in <ref>.

§.§ Parity of Qualified Persons

The following probabilistic verification problem concerns verifying the parity of qualified persons, a variant of demographic parity that only considers the subpopulation of persons qualified for, for example, hiring <cit.>.
Let \(\mathcal{X} \subset \Reals^n\), \(A \subset \{1, \ldots, n\}\), \(a \in A\), and \(\NN: \Reals^n \to \Reals^2\) be as in <ref>.
Additionally, let \(q \in \{1, \ldots, n\} \setminus A\) and \(\hat{q} \in \Reals\), such that persons with \(\inp_q \geq \hat{q}\) are considered to be qualified.
In their extended set of experiments, AlbarghouthiDAntoniDrewsEtAl2017 consider a \(q\) that encodes age and \(\hat{q} = 18\) so that only persons who are at least 18 years old are considered to be qualified.
The parity of qualified persons fairness notion is

[x⃗](() =  | is disadvantaged ∧ is qualified)
[x⃗](() =  | is advantaged ∧ is qualified)
≥    γ

[](()_1 - ()_2 ≥0 |_a ≤0 ∧_q ≥q̂)
[](()_1 - ()_2 ≥0 |_a ≥1 ∧_q ≥q̂)
≥    γ

[](()_1 - ()_2 ≥0 |min(-_a, _q - q̂) ≥0)
[](()_1 - ()_2 ≥0 |min(_a - 1, _q - q̂) ≥0)
≥    γ

⟺    (
[](*[1](, ()) ≥0),
[](*[4](, ()) ≥0)
)    ≥    0

where \(\gamma \in [0, 1]\), \(\fSAT(p_1,p_2,p_3,p_4) = (p_1 p_4)/(p_2 p_3) - \gamma\), \(\fSAT*[1](\inp, \NN(\inp)) = \min({\NN(\inp)}_1 - {\NN(\inp)}_2, -\inp_a, \inp_q - \hat{q})\), \(\fSAT*[2](\inp, \NN(\inp)) = \min(-\inp_a, \inp_q - \hat{q})\), \(\fSAT*[3](\inp, \NN(\inp)) = \min({\NN(\inp)}_1 - {\NN(\inp)}_2, \inp_a - 1, \inp_q - \hat{q})\), and \(\fSAT*[4](\inp, \NN(\inp)) = \min(\inp_a -1, \inp_q - \hat{q})\).

§.§ ACAS Xu Safety

Next, we consider <ref> for an ACAS Xu network, where to be safe means satisfying property \(\phi_2\) of KatzBarrettDillEtAl2017.
For quantifying the number of violations, we first define what it means for an ACAS Xu neural network \(\NN: \Reals^5 \to \Reals^5\) to violate \(\phi_2\).
Using the satisfaction functions of BauerMarquartBoetiusLeueEtAl2021, violating \(\phi_2\) means

\begin{equation}
\fSAT*(\inp, \NN(\inp)) = \max_{i = 1}^{5} {\NN(\inp)}_i - {\NN(\inp)}_1 < 0
\quad \forall \inp \in \mathcal{X}_{\phi_2} \cap \mathcal{X},\label{eqn:acasxu-phitwo-verif}
\end{equation}

where \(\mathcal{X}\) is the bounded hyperrectangular input space of \(\NN\) and

\begin{equation*}
\mathcal{X}_{\phi_2} = [55947.961, \infty] \times \Reals^2 \times [1145, \infty] \times [-\infty, 60].
\end{equation*}

We refer to KatzBarrettDillEtAl2017 for an interpretation of \(\phi_2\) in the application context.
Quantifying the number of violations with respect to \(\phi_2\) corresponds to computing

\begin{equation*}
\ell \leq \prob[\inp](\fSAT*(\inp, \NN(\inp)) < 0) = \prob[\inp](-\fSAT*(\inp, \NN(\inp)) \geq 0) \leq u,
\end{equation*}

where \(\fSAT*\) is as in <ref> and \(\inp\) is uniformly distributed on \(\mathcal{X}_{\phi_2} \cap \mathcal{X}\) with all points outside \(\mathcal{X}_{\phi_2} \cap \mathcal{X}\) having zero probability.

§.§ ACAS Xu Robustness

For the ACAS Xu robustness experiment in <ref>, we solve five probabilistic verification problems for each reference input \(\inp\) — one for each of the five classes.
Our goal is to bound the probability of \(\NN\) classifying an input \(\vec{x}'\) as class \(i \in \{1, \ldots, 5\}\), where \(\vec{x}'\) is close to the reference input \(\inp\) in the first two dimensions and identical to \(\inp\) in the remaining dimensions.

Let \(\NN\) be the ACAS Xu network \(N_{1,1}\) of KatzBarrettDillEtAl2017 with input space \(\mathcal{X} = \HR{\inp}\).
Let \(\inp\) be a reference input.
Note that the ACAS Xu networks assign the class with the minimal score to an input instead of using the maximal score.
For bounding the probability of obtaining class \(i \in \{1, \ldots, 5\}\) for inputs close to \(\inp\), we compute bounds on

\begin{equation*}
\prob[\vec{x}'](\fSAT*(\vec{x}', \NN(\vec{x}')) \geq 0),
\end{equation*}

where \(\fSAT*(\vec{x}', \NN(\vec{x}')) = \min_{j = 1}^5 {\NN(\vec{x}')}_j - {\NN(\vec{x}')}_i\) and \(\vec{x}'\) is uniformly distributed on the set

\begin{equation*}
\mathcal{X} \cap
\left(
[\inp_{1:2} - 0.05\cdot\vec{w}_{1:2}, \inp_{1:2} + 0.05\cdot\vec{w}_{1:2}]
\times \{\inp_{3:5}\}
\right),
\end{equation*}

where \(\vec{w} = \overline{\inp} - \underline{\inp}\) and \(\vec{z}_{i:j}\) is the vector containing the elements \(i, \ldots, j\) of a vector \(\vec{z}\).

§ INTERVAL ARITHMETIC

This section introduces additional interval arithmetic bounding rules for linear functions, multiplication, and division, complementing the interval arithmetic bounding rules for monotone functions in <ref>.
Furthermore, we provide Theorems 5.1 and 6.1 of MooreKearfottCloud2009 for reference.

§.§ Bounding Rules

Let \(f^{(k)}\) be as in <ref>.
First, consider the multiplication of two scalars, that is, \(f^{(k)}(z, w) = zw\) where \(\underline{z} \leq z \leq \overline{z}\) and \(\underline{w} \leq w \leq \overline{w}\).
We have

\begin{equation*}
\min(\underline{z}\underline{w}, \underline{z} \overline{w}, \overline{z}\underline{w}, \overline{z}\overline{w})
\leq z_1z_2 \leq
\max(\underline{z}\underline{w}, \underline{z} \overline{w}, \overline{z}\underline{w}, \overline{z}\overline{w}).
\end{equation*}

For the element-wise multiplication of vectors, we apply the above rule to each element separately.
Multiplication of several arguments can be rewritten as several multiplications of two arguments.

Now, consider computing bounds of the reciprocal \(f^{(k)}(z) = \frac{1}{z}\) with \(\underline{z} \leq z \leq \overline{z}\).
We differentiate the following cases
1/z ≤1/z     if  0 ∉(z, z]

1/z ≤1/z     if  0 ∉[z, z)

-∞≤1/z     if  0 ∈(z, z]

1/z ≤∞    if  0 ∈[z, z).

Using bounds on the reciprocal, we can compute bounds on a division by rewriting division as multiplication by the reciprocal.
Lastly, for a linear function \(f^{(k)}(\vec{z}) = \mat{A}\vec{z} + \vec{b}\) where \(\underline{\vec{z}} \leq \vec{z} \leq \overline{\vec{z}}\), we have

\begin{equation*}
\mat{A}^+\underline{z} + \mat{A}^-\overline{z} + \vec{b}
\leq \mat{A} \vec{z} + \vec{b} \leq
\mat{A}^+\overline{z} + \mat{A}^-\underline{z} + \vec{b},
\end{equation*}

where \(\mat{A}^+_{i,j} = \max(0, \mat{A}_{i,j})\) and \(\mat{A}^-_{i,j} = \min(0, \mat{A}_{i,j})\).

§.§ Selected Theorems of MooreKearfottCloud2009

We include Theorems 5.1 and 6.1 of MooreKearfottCloud2009 and relevant definitions for reference.
Let \(\mathbb{H}^n = \{[\underline{\vec{x}}, \overline{\vec{x}}] \mid \underline{\vec{x}}, \overline{\vec{x}} \in \Reals^n, \underline{\vec{x}} \leq \overline{\vec{x}}\}\) be the set of hyperrectangles in \(\Reals^n\).
Let \(w: \powerset{\Reals^n} \to \RealsNonNeg\) with

\begin{equation*}
w(\mathcal{X}') = \max_{i \in \{1, \ldots n\}}{(\max_{\vec{x} \in \mathcal{X}'}\vec{x}_i - \min_{\vec{x} \in \mathcal{X}'}\vec{x}_i)}
\end{equation*}

be the width of a set.
We denote the image of a hyperrectangle \(\HR{\inp}\) under \(f: \Reals^n \to \Reals^m\) as \(f(\HR{\inp}) = \{f(\inp) \mid \inp \in \HR{\inp}\}\).

§.§.§ Definitions

Theorem 5.1 of MooreKearfottCloud2009 applies to inclusion isotonic interval extensions as defined below.

Let \(F: \mathbb{H}^n \to \mathbb{H}^m\) and \(f: \Reals^n \to \Reals^m\).

* \(F\) is an interval extensions of \(f\) if \(\forall \vec{x} \in \Reals^n: F([\vec{x}, \vec{x}]) = [f(\vec{x}), f(\vec{x})]\).

* \(F\) is inclusion isotonic if \(\forall [\underline{\vec{x}}, \overline{\vec{x}}], [\underline{\vec{x}}', \overline{\vec{x}}'] \in \mathbb{H}^n, [\underline{\vec{x}}', \overline{\vec{x}}'] \subseteq [\underline{\vec{x}}, \overline{\vec{x}}]: F([\underline{\vec{x}}', \overline{\vec{x}}']) \subseteq F([\underline{\vec{x}}, \overline{\vec{x}}])\).

Theorem 6.1 of MooreKearfottCloud2009 requires an inclusion isotonic Lipschitz interval extension.

Let \(F: \mathbb{H}^n \to \mathbb{H}^m\) be an interval extension of \(f: \Reals^n \to \Reals^m\).
\(F\) is Lipschitz if there exists an \(L \in \Reals\) such that \(\forall [\underline{\vec{x}}, \overline{\vec{x}}] \in \mathbb{H}^n: w(F([\underline{\vec{x}}, \overline{\vec{x}}])) \leq Lw(\HR{\vec{x}})\).

Interval arithmetic, as introduced in <ref>, corresponds to natural interval extensions in MooreKearfottCloud2009.
As MooreKearfottCloud2009 show, natural interval extensions — and, therefore, interval arithmetic as introduced in <ref> — satisfy Definitions <ref> and <ref>.

§.§.§ Theorems

Theorem 5.1 of MooreKearfottCloud2009 is known as the fundamental theorem of interval analysis.
Theorem 6.1 of MooreKearfottCloud2009 underlines our completeness proof in <ref>.

[Theorem 5.1 of MooreKearfottCloud2009]
If \(F: \mathbb{H}^n \to \mathbb{H}^m\) is an inclusion isotonic interval extension of \(f: \Reals^n \to \Reals^m\), we have \(f(\HR{\inp}) \subseteq F(\HR{\inp})\) for every \(\HR{\inp} \in \mathbb{H}^n\).

[Theorem 6.1 of MooreKearfottCloud2009]
Let \(F: \mathbb{H}^n \to \mathbb{H}^m\) be an inclusion isotonic Lipschitz interval extension of \(f: \Reals^n \to \Reals^m\).
Let \(\mathcal{X} = \HR{\inp} \in \mathbb{H}^n\).
We define the \(M\)-step uniform subdivision of \(\HR{\inp}\) with \(M \in \Nats\) as

\begin{equation*}
\mathcal{X}_{i,j} = \left[
\underline{\inp}_i + (j - 1)\frac{w([\underline{\inp}_i, \overline{\inp}_i])}{M},
\underline{\inp}_i + j\frac{w([\underline{\inp}_i, \overline{\inp}_i])}{M}
\right],
\quad j \in \{1, \ldots, M\}.
\end{equation*}

Further, let

\begin{equation*}
F^{(M)}(\HR{\inp}) = \bigcup_{j_i = 1}^{M} F(\mathcal{X}_{1,j_1} \times \cdots \times \mathcal{X}_{n,j_n}). \\
\end{equation*}

It holds that

\begin{equation*}
w(F^{(M)}(\HR{\inp})) - w(f(\HR{\inp})) \leq 2L\frac{w(\mathcal{X})}{M},
\end{equation*}

where \(L\) is the Lipschitz constant of \(F\).

In this paper, we apply Theorem 6.1 to a refinement of the input space finer than some uniform subdivision, meaning that every part of the input space refinement is contained in a part of the uniform subdivision.
Theorem 6.1 extends to such refinements of the input space since interval arithmetic is inclusion isotonic.

§ EXPERIMENTS

This section contains additional details on the experiments from <ref> and additional experimental results.

§.§ Hardware

We run all our experiments in <ref> and this section on a workstation running Ubuntu 22.04 with an Intel i7–4820K CPU, 32 GB of memory and no GPU (HW1).
This CPU model is ten years old (introduced in late 2013) and has four cores and eight threads.
ConverseFilieriGopinathEtAl2020 use an AMD EPYC 7401P CPU for their ACAS Xu robustness experiments, limiting their tool to use 46 threads and at most 4GB of memory.
AMD EPYC 7401P was introduced in mid-2017, targeting servers.
MarzariRoncolatoFarinelli2023 use a Ubuntu 22.04 workstation with an Intel i5–13600KF CPU and an Nvidia GeForce RTX 4070 Ti GPU.
This CPU was introduced in end-2022 and has 14 cores with 20 threads.
<Ref> contains the CPU comparison from <versus.com>.
As the figure shows, our HW1 is less performant when considering both the theoretical performance according to the hardware specification and the actual performance on a series of computational benchmarks.
Therefore, HW1 is inferior in computation power to the hardware used by ConverseFilieriGopinathEtAl2020 and MarzariRoncolatoFarinelli2023.
Additional experiments on recent hardware are contained in <ref>.

HW1 VersusComOurs vs. ConverseFilieriGopinathEtAl2020 VersusComTheirs

HW1 VersusComOurs vs. MarzariRoncolatoFarinelli2023 VersusComTheirs

Hardware Comparison from <versus.com>.
The figures are taken from <https://versus.com/en/amd-epyc-7401p-vs-intel-core-i7-6700>
and <https://versus.com/en/intel-core-i5-13600kf-vs-intel-core-i7-4820k>, respectively.
Printouts of both web pages are available in the supplementary material.

§.§ FairSquare Benchmark

We provide additional details on the FairSquare benchmark and present additional results.

§.§.§ Extended Description of the Benchmark

The input space in the FairSquare benchmark consists of three continuous variables for the age, the years of education (edu), and the yearly gain in capital of a person (capital gain), respectively, as well as an additional discrete protected variable indicating a person's (assumed binary) sex.
The neural networks use age and edu or age, edu and capital gain, depending on the network architecture, to predict whether a person has a high salary (higher than $50 000).
The three input distributions used by the FairSquare benchmark are

* the combination of three independent normal distributions for the continuous variables and a Bernoulli distribution for the person's sex,

* a Bayesian Network with the structure in <ref> introducing correlations between the variables that are similarly distributed as for the independent input distribution (Bayes Net 1), and

* the same Bayesian Network augmented with a constraint that the years of education may not exceed a person's age (Bayes Net 2).

The integrity constraint is implemented by AlbarghouthiDAntoniDrewsEtAl2017 as a post-processing of the samples from the Bayesian Network.
In particular, an person's years of education are set to the person's age if the sampled years of education are larger than the sampled age.
We introduce a pre-processing layer before the neural network we want to verify to implement this integrity constraint in our setting.
This pre-processing layer computes \(\mathrm{edu}' = \min(\mathrm{edu}, \mathrm{age})\) and leaves the remaining inputs unchanged.

Besides experiments on the demographic parity fairness notion, the extended FairSquare benchmarks also include experiments on the parity of qualified persons fairness notion as defined in <ref>.
For both fairness notions, the FairSquare benchmark uses a fairness threshold of \(\gamma = 0.85\).

[minimum size=1cm]
[draw=black, circle] (S) at (0,0) sex;
[draw=black, circle, below=of S] (C) cg;
[draw=black, circle, below left=of C] (A) age;
[draw=black, circle, below right=of C] (E) edu;

[black,-Stealth] (S) – (C);
[black,-Stealth] (C) – (A);
[black,-Stealth] (C) – (E);

FairSquare Bayes Net 1.
The network structure of the Bayesian Network from the FairSquare benchmark.
In this figure, cg denotes the capital gain variable and edu denotes the number of years of education.

§.§.§ Extended Results

<Ref> contains the comparison of  and  on the extended set of FairSquare benchmarks from AlbarghouthiDAntoniDrewsEtAl2017.
The extended set of benchmarks also contains one instance which  can not solve within the time budget.
In this instance, the bounds computed up to the timeout indicate that the probability ratio on which we compute bounds (see <ref>) is close to the threshold of \(0.85\) above which the probability ratio should lie according to the fairness notion.
Relating this to <ref>, the above might mean that the neural network is not strictly feasible or infeasible for the probabilistic verification problem, at least not by a large margin.
This shows that the scalability of  not only depends on the input size and the network complexity but also on additional factors, such as the margin with which a neural network is feasible or infeasible for the probabilistic verification problem.

As <ref> shows, there are also two benchmark instances in which  is substantially faster than .
However, these are outliers in the entire set of benchmarks.
Additionally, they disappear when considering input distributions that were clipped to realistic ranges for the inputs (for example, no negative values of age).
<Ref> contains the results of running  on these clipped input distributions.
Using realistic input distributions accelerates , except for three benchmark instances.
It also reveals that some FairSquare networks are, in fact, unfair under realistic conditions.
Faithfully comparing these results with  is not feasible since  does not implement truncated normal distributions, which we use for the clipped input distributions.

3cParity of

3cDemographic Parity
3cQualified Persons

2cBenchmark Instance
2cRuntime (s)
2cRuntime (s)

$\NN$     $\prob[\inp]$
FairSquare     Fair?
FairSquare     Fair?
NN2,1     independent     1.3     2.1            1.2     3.2

NN2,1     Bayes Net 1     2.2     75.5           3.2     37.8

NN2,1     Bayes Net 2     3.6     50.5           6.1     25.5

NN2,2     independent     1.8     4.2            1.7     6.8

NN2,2     Bayes Net 1     3.3     33.3           622.2            60.4

NN2,2     Bayes Net 2     4.9     58.6           782.8            60.7

NN3,2     independent     2.2     450.5          2.0     657.0

NN3,2     Bayes Net 1     4.4     TO             TO             TO

NN3,2     Bayes Net 2     6.9     TO             11.5      TO

FairSquare Benchmark.
The first two columns indicate the neural network \(\NN\) that is verified and the input probability distribution \(\prob[\inp]\), respectively.
We verify two fairness notions: demographic parity and parity of qualified persons.
For each fairness notion, the first two columns contain the runtime of  and  in seconds.
Here, TO indices timeout (900s).
The last column contains the result of the verification with .
The results in this table are for the non-clipped input distributions that include unrealistic values, such as negative values of a person's age.

2cParity of

2cDemographic Parity
2cQualified Persons

2cBenchmark Instance

$\NN$     $\prob[\inp]$
Runtime     Fair?
Runtime     Fair?
NN2,1     independent clipped     1.1          1.1

NN2,1     Bayes Net 1 clipped     1.8          2.6

NN2,1     Bayes Net 2 clipped     2.9          5.4

NN2,2     independent clipped     1.7          1.6

NN2,2     Bayes Net 1 clipped     2.8          3.8

NN2,2     Bayes Net 2 clipped     4.3          7.5

NN3,2     independent clipped     1.3          1.3

NN3,2     Bayes Net 1 clipped     5.0          6.2

NN3,2     Bayes Net 2 clipped     7.8           TO

FairSquare Benchmark for Realistic Input Distributions.
This table contains the results for running  on the FairSquare benchmarks with clipped input distributions that rule out unrealistic values, such as negative values of a person's age.
The first two columns indicate the neural network \(\NN\) that is verified and the input probability distribution \(\prob[\inp]\), respectively.
The following columns contain the runtime of  in seconds and the verification result for the demographic parity and the parity of qualified persons fairness notions.
A runtime of TO indices timeout (900s).

§.§ ACAS Xu Safety

We provide a comparison of  with  <cit.> on all 36 ACAS Xu networks to which property \(\phi_2\) applies <cit.>.
These results are contained in <ref>, revealing that the networks in <ref> are outliers.
However,  still computes tighter sound bounds faster than  computes probably sound bounds for 13 of the 36 networks.

<Ref> contains additional results for running  with a finer grid of time budgets on the networks from <ref> and two additional challenging instances from KatzBarrettDillEtAl2017.

6c     3c
2c10s     2c30s     2c60s     3c99.9% confid.

$\phi$     $\NN$     $\ell,u$         $\ell,u$         $\ell,u$         $u$         1cRt

$N_{2,1}$     $  0.04\%,   3.14\%$             $  0.08\%,   2.48\%$    $\ms$    $  0.16\%,   2.05\%$            $  2.58\%$             64.00s

$N_{2,2}$     $  0.00\%,   3.55\%$             $  0.11\%,   2.92\%$             $  0.24\%,   2.67\%$            $  2.49\%$    $\ms$     43.60s

$N_{2,3}$     $  0.54\%,   3.43\%$    $\ms$    $  0.89\%,   3.05\%$             $  1.01\%,   2.69\%$            $  3.58\%$             54.50s

$N_{2,4}$     $  0.00\%,   2.74\%$    $\ms$    $  0.03\%,   2.30\%$             $  0.05\%,   2.08\%$            $  2.78\%$             65.90s

$N_{2,5}$     $  0.05\%,   4.16\%$             $  0.19\%,   3.70\%$             $  0.41\%,   3.22\%$            $  3.06\%$    $\ms$     47.20s

$N_{2,6}$     $  0.02\%,   3.71\%$             $  0.20\%,   3.09\%$             $  0.33\%,   2.58\%$            $  2.32\%$    $\ms$     45.40s

$N_{2,7}$     $  0.42\%,   5.54\%$             $  1.02\%,   5.03\%$             $  1.32\%,   4.42\%$            $  4.04\%$    $\ms$     47.40s

$N_{2,8}$     $  0.98\%,   3.80\%$             $  1.25\%,   3.04\%$    $\ms$    $  1.38\%,   2.79\%$            $  3.28\%$             52.90s

$N_{2,9}$     $  0.02\%,   2.90\%$             $  0.09\%,   1.67\%$             $  0.11\%,   1.41\%$            $  0.58\%$    $\ms$     10.80s

$N_{3,1}$     $  0.26\%,   4.26\%$             $  0.55\%,   3.33\%$             $  0.88\%,   2.78\%$    $\ms$    $  2.84\%$             40.30s

$N_{3,2}$     $  0.00\%,   3.45\%$             $  0.00\%,   1.58\%$             $  0.00\%,   1.16\%$            $  0.00\%$    $\ms$      1.00s

$N_{3,3}$     $  0.00\%,   3.06\%$             $  0.00\%,   1.84\%$             $  0.00\%,   0.87\%$            $  0.00\%$    $\ms$      1.00s

$N_{3,4}$     $  0.00\%,   3.24\%$             $  0.07\%,   1.96\%$             $  0.13\%,   1.52\%$            $  1.36\%$    $\ms$     30.80s

$N_{3,5}$     $  0.37\%,   2.83\%$             $  0.54\%,   2.30\%$    $\ms$    $  0.66\%,   2.12\%$            $  2.67\%$             41.90s

$N_{3,6}$     $  0.00\%,   6.73\%$             $  0.00\%,   5.50\%$             $  0.03\%,   5.26\%$            $  2.69\%$    $\ms$     31.60s

$N_{3,7}$     $  0.00\%,   5.77\%$             $  0.00\%,   4.61\%$             $  0.00\%,   3.66\%$            $  0.91\%$    $\ms$     30.10s

$N_{3,8}$     $  0.17\%,   5.48\%$             $  0.27\%,   4.68\%$             $  0.34\%,   3.76\%$            $  2.79\%$    $\ms$     61.40s

$N_{3,9}$     $  0.34\%,   4.77\%$             $  1.07\%,   4.02\%$             $  1.42\%,   3.73\%$            $  3.52\%$    $\ms$     32.10s

$N_{4,1}$     $  0.00\%,   2.89\%$             $  0.02\%,   1.95\%$             $  0.05\%,   1.57\%$            $  1.16\%$    $\ms$     26.40s

$N_{4,2}$     $  0.00\%,   2.55\%$             $  0.00\%,   1.65\%$             $  0.00\%,   1.33\%$            $  0.00\%$    $\ms$      1.00s

$N_{4,3}$     $  0.19\%,   2.90\%$             $  0.37\%,   2.66\%$             $  0.62\%,   2.28\%$    $\ms$    $  2.43\%$             40.80s

$N_{4,4}$     $  0.05\%,   2.42\%$             $  0.18\%,   2.23\%$             $  0.28\%,   2.00\%$            $  1.94\%$    $\ms$     40.00s

$N_{4,5}$     $  0.04\%,   4.66\%$             $  0.23\%,   3.58\%$             $  0.64\%,   3.02\%$            $  2.85\%$    $\ms$     40.70s

$N_{4,6}$     $  0.39\%,   4.45\%$             $  0.97\%,   3.52\%$             $  1.22\%,   3.28\%$            $  3.22\%$    $\ms$     38.80s

$N_{4,7}$     $  0.18\%,   4.28\%$             $  0.39\%,   3.42\%$             $  0.63\%,   3.02\%$            $  2.58\%$    $\ms$     29.60s

$N_{4,8}$     $  0.20\%,   3.89\%$             $  0.53\%,   3.27\%$    $\ms$    $  0.88\%,   2.82\%$            $  3.51\%$             52.40s

$N_{4,9}$     $  0.00\%,   3.39\%$             $  0.00\%,   2.08\%$             $  0.01\%,   1.65\%$            $  0.39\%$    $\ms$    11.00s

$N_{5,1}$     $  0.03\%,   2.57\%$             $  0.32\%,   2.13\%$    $\ms$    $  0.44\%,   2.00\%$            $  2.14\%$            34.80s

$N_{5,2}$     $  0.24\%,   2.71\%$    $\ms$    $  0.49\%,   2.06\%$             $  0.53\%,   1.95\%$            $  3.22\%$            66.80s

$N_{5,3}$     $  0.00\%,   1.74\%$             $  0.00\%,   0.73\%$             $  0.00\%,   0.53\%$            $  0.00\%$    $\ms$     1.00s

$N_{5,4}$     $  0.01\%,   2.95\%$             $  0.15\%,   2.35\%$             $  0.24\%,   2.09\%$    $\ms$    $  2.31\%$            54.20s

$N_{5,5}$     $  0.77\%,   3.37\%$             $  1.33\%,   2.68\%$    $\ms$    $  1.38\%,   2.63\%$            $  2.77\%$            32.90s

$N_{5,6}$     $  0.34\%,   4.89\%$             $  0.87\%,   3.36\%$    $\ms$    $  1.05\%,   3.07\%$            $  4.82\%$            73.50s

$N_{5,7}$     $  1.29\%,   4.62\%$             $  1.82\%,   4.12\%$    $\ms$    $  2.00\%,   3.75\%$            $  4.27\%$            44.60s

$N_{5,8}$     $  0.95\%,   4.12\%$             $  1.37\%,   3.33\%$    $\ms$    $  1.58\%,   3.11\%$            $  3.38\%$            42.40s

$N_{5,9}$     $  0.78\%,   3.74\%$             $  1.28\%,   3.06\%$    $\ms$    $  1.46\%,   2.91\%$            $  3.06\%$            35.60s

2r$\#\!\ms$     2r3     2r10     2r3     2r20

[Comparison of ProbBounds and eProVe for ACAS Xu Safety]
Comparison of  and  for ACAS Xu Safety.
A  column is marked with an arrow (\(\ms\)) if  computes a tighter upper bound than  within a certain time budget.
If this is not the case for any time budget, the  entry is marked with an arrow (\(\ms\)).
In total,  computes tighter upper bounds within 60s in 16 cases.
In 13 cases, it even computes a tighter bound faster than .

2c10s     2c30s     2c1m
(lr)3-4 (lr)5-6 (lr)7-8
$\phi$     $\NN$     $\ell,u$     $u-\ell$     $\ell,u$     $u-\ell$     $\ell,u$     $u-\ell$

$N_{4,3}$     $  0.19\%,   2.90\%$     $  2.71\%$     $  0.37\%,   2.66\%$     $  2.29\%$     $  0.62\%,   2.28\%$     $  1.66\%$

$N_{4,9}$     $  0.00\%,   3.39\%$     $  3.39\%$     $  0.00\%,   2.08\%$     $  2.08\%$     $  0.01\%,   1.65\%$     $  1.64\%$

$N_{5,8}$     $  0.95\%,   4.12\%$     $  3.17\%$     $  1.37\%,   3.33\%$     $  1.97\%$     $  1.58\%,   3.11\%$     $  1.53\%$

$N_{1,9}$     $  0.00\%,  94.83\%$     $ 94.83\%$     $  0.00\%,  86.28\%$     $ 86.28\%$     $  0.00\%,  77.96\%$     $ 77.96\%$

$N_{2,9}$     $  0.00\%,  70.06\%$     $ 70.06\%$     $  0.00\%,  59.38\%$     $ 59.38\%$     $  0.00\%,  53.51\%$     $ 53.51\%$

2c10m     2c1h
$\phi$     $\NN$     $\ell,u$     $u-\ell$     $\ell,u$     $u-\ell$

$N_{4,3}$     $  0.97\%,   1.97\%$     $  1.01\%$     $  1.12\%,   1.75\%$     $  0.63\%$

$N_{4,9}$     $  0.04\%,   0.55\%$     $  0.51\%$     $  0.08\%,   0.30\%$     $  0.22\%$

$N_{5,8}$     $  1.91\%,   2.66\%$     $  0.76\%$     $  1.97\%,   2.57\%$     $  0.60\%$

$N_{1,9}$     $  0.00\%,  44.48\%$     $ 44.48\%$     $  0.00\%,  27.66\%$     $ 27.66\%$

$N_{2,9}$     $  0.00\%,  31.57\%$     $ 31.57\%$     $  0.00\%,  14.73\%$     $ 14.73\%$

Extended  Results for ACAS Xu Safety.

§.§ ACAS Xu Robustness

<Ref> contains the bounds computed by  for each reference input and each output class: COC (Clear-of-Conflict), WL (steer Weak Left), WR (steer Weak Right), SL (steer Strong Left), and SR (steer Strong Right).
The table reveals that the ACAS Xu network \(N_{1,1}\) tends to classify inputs as Clear-of-Conflict (COC), regardless of the class assigned to the reference input.
This insight does not agree with the insight drawn by ConverseFilieriGopinathEtAl2020.
However, both results are based on a tiny sample of the input space of only 25 points.
Obtaining valid results requires considering a significantly larger sample of the input space or quantifying global robustness <cit.>.

10cTarget Label
2cReference Input     2cCOC     2cWL     2cWR     2cSL     2cSR
Label     Input #
$\ell, u$     1cRt     $\ell, u$     1cRt     $\ell, u$     1cRt     $\ell, u$     1cRt     $\ell, u$     1cRt
COC     1     $ 100\%,  100\%$        0.3     $0.0\%, 0.0\%$        0.3     $0.0\%, 0.0\%$        0.3     $0.0\%, 0.0\%$        0.3     $0.0\%, 0.0\%$        0.3

COC     2     $99.9\%,  100\%$        0.6     $0.0\%, 0.0\%$        0.5     $0.0\%, 0.1\%$        0.5     $0.0\%, 0.0\%$        0.5     $0.0\%, 0.0\%$        0.5

COC     3     $91.7\%, 91.8\%$        2.9     $1.3\%, 1.4\%$        5.4     $1.9\%, 2.0\%$        3.9     $2.1\%, 2.2\%$       10.1     $2.8\%, 2.8\%$        3.4

COC     4     $ 100\%,  100\%$        0.2     $0.0\%, 0.0\%$        0.2     $0.0\%, 0.0\%$        0.2     $0.0\%, 0.0\%$        0.2     $0.0\%, 0.0\%$        0.2

COC     5     $ 100\%,  100\%$        0.4     $0.0\%, 0.0\%$        0.4     $0.0\%, 0.0\%$        0.4     $0.0\%, 0.0\%$        0.4     $0.0\%, 0.0\%$        0.4

WL      1     $89.9\%, 90.0\%$        1.6     $1.6\%, 1.7\%$       30.8     $4.7\%, 4.8\%$        1.8     $3.6\%, 3.7\%$       31.8     $0.0\%, 0.0\%$        1.4

WL      2     $92.9\%, 93.0\%$        2.9     $2.4\%, 2.4\%$        2.0     $3.0\%, 3.0\%$        2.7     $1.4\%, 1.5\%$        2.5     $0.2\%, 0.3\%$        1.4

WL      3     $90.4\%, 90.5\%$        2.9     $1.1\%, 1.2\%$        7.5     $3.6\%, 3.7\%$        3.3     $3.8\%, 3.9\%$        9.1     $0.8\%, 0.9\%$        2.1

WL      4     $83.6\%, 83.7\%$       16.5     $5.4\%, 5.5\%$       44.2     $2.6\%, 2.7\%$       41.1     $4.6\%, 4.7\%$       49.3     $3.5\%, 3.6\%$       51.2

WL      5     $96.8\%, 96.9\%$        1.5     $2.9\%, 3.0\%$        1.4     $0.0\%, 0.1\%$        1.0     $0.0\%, 0.0\%$        0.8     $0.1\%, 0.2\%$        1.0

WR      1     $62.2\%, 62.3\%$       15.2     $0.0\%, 0.0\%$        1.1     $3.4\%, 3.5\%$       16.1     $15.8\%, 15.9\%$        7.3     $18.5\%, 18.6\%$        9.9

WR      2     $96.6\%, 96.7\%$        1.3     $1.5\%, 1.6\%$        1.3     $1.7\%, 1.8\%$        1.5     $0.0\%, 0.1\%$        1.3     $0.0\%, 0.1\%$        1.0

WR      3     $94.8\%, 94.9\%$        3.1     $0.0\%, 0.1\%$       24.1     $1.9\%, 2.0\%$      173.4     $0.0\%, 0.1\%$        4.7     $3.2\%, 3.3\%$       92.9

WR      4     $94.3\%, 94.4\%$        1.8     $1.8\%, 1.9\%$        2.2     $3.4\%, 3.5\%$        1.9     $0.3\%, 0.4\%$        2.0     $0.0\%, 0.1\%$        1.2

WR      5     $91.1\%, 91.2\%$        3.8     $0.8\%, 0.9\%$        6.3     $4.2\%, 4.3\%$        7.5     $3.0\%, 3.1\%$       10.0     $0.7\%, 0.8\%$        4.1

SL      1     $81.1\%, 81.2\%$       20.2     $2.3\%, 2.4\%$       60.1     $4.8\%, 4.9\%$      117.9     $4.8\%, 4.9\%$       58.3     $6.7\%, 6.8\%$      118.5

SL      2     $93.9\%, 93.9\%$        2.7     $1.1\%, 1.2\%$       21.2     $3.7\%, 3.8\%$       31.3     $0.7\%, 0.7\%$        9.2     $0.5\%, 0.6\%$       12.4

SL      3     $83.3\%, 83.4\%$        9.1     $2.8\%, 2.9\%$       13.7     $6.5\%, 6.6\%$       33.5     $2.1\%, 2.2\%$       12.7     $5.0\%, 5.1\%$       27.9

SL      4     $81.8\%, 81.9\%$       18.0     $3.3\%, 3.4\%$       80.0     $5.2\%, 5.3\%$       99.0     $5.6\%, 5.7\%$       77.5     $4.0\%, 4.1\%$       58.8

SL      5     $84.9\%, 85.0\%$       13.3     $1.7\%, 1.8\%$       22.4     $5.4\%, 5.5\%$       35.5     $3.0\%, 3.1\%$       23.7     $4.8\%, 4.9\%$       34.3

SR      1     $91.0\%, 91.1\%$        4.9     $0.0\%, 0.1\%$        2.6     $5.2\%, 5.3\%$       50.8     $0.0\%, 0.1\%$        2.2     $3.6\%, 3.7\%$       28.5

SR      2     $87.0\%, 87.1\%$        7.3     $0.9\%, 1.0\%$        5.0     $5.7\%, 5.8\%$       56.4     $1.2\%, 1.3\%$        6.1     $5.0\%, 5.1\%$       72.3

SR      3     $88.4\%, 88.5\%$        2.7     $2.5\%, 2.6\%$        8.1     $4.7\%, 4.8\%$        7.3     $2.6\%, 2.7\%$       12.1     $1.5\%, 1.6\%$        7.4

SR      4     $93.0\%, 93.1\%$        5.4     $0.7\%, 0.8\%$       47.6     $3.7\%, 3.8\%$       86.9     $1.0\%, 1.1\%$       25.5     $1.4\%, 1.5\%$       29.6

SR      5     $79.3\%, 79.4\%$        5.2     $0.2\%, 0.3\%$        1.7     $9.1\%, 9.2\%$        8.3     $3.7\%, 3.8\%$        4.2     $7.4\%, 7.5\%$        5.1

ACAS Xu Robustness.
We report the lower and upper bounds (\(\ell, u\))  computes for the probability that a perturbed input is classified as the target label and the runtime in seconds (Rt) of  for computing these bounds.
For each combination of unperturbed label, input, and target label, we run  until the difference between the lower and the upper bound is at most \(0.1\%\).

§.§ MiniACSIncome

We provide additional details on the MiniACSIncome benchmark, including how we construct the dataset and train the networks.

§.§.§ Dataset

The MiniACSIncome benchmarks are built by sampling about 100 000 entries from the ACS PUMPS 1-Year horizon data for all states of the USA for the year 2018 using the folktables Python package <cit.>.
In line with ACSIncome <cit.>, we only sample individuals older than 16 years, with a yearly income of at least $100, reported working hours per week of at least 1, and a PWGTP (more details in DingHardtMillerEtAl2021) of at least 1.
In total, our dataset contains 102 621 samples.

Variable Order.

For obtaining the benchmark MiniACSIncome-$i$, \(i \in \{1, \ldots, 8\}\), we select \(i\) input variables in the following order: SEX, COW, SCHL, WKHP, MAR, RAC1P, RELP, AGEP.
We choose SEX as the first variable so that we can verify the fairness with respect to SEX on every benchmark instance.
The order of the remaining variables is chosen based on each variable's expected predictive value and the number of discrete values.
In particular, we select COW (class of work), SCHL (level of education), and WKHP (work hours per week) first, as we consider these variables to be more predictive than MAR (marital status), RAC1P (races of a person), RELP (relationship), and AGEP (age of a person).
The variables are ordered by their number of discrete values within these groups of expected predictive value.
For example, COW has nine categories, while WKHP has 99 possible integer values.

<Ref> contains the number of input dimensions and the total number of discrete values in each MiniACSIncome-$i$ input space.
The input space contains more dimensions than input variables due to the one-hot encoding of all categorical variables.

8c#Input Variables

1     2     3     4     5     6     7     8
#Input Dimensions     2     10     34     35     40     49     67     68

#Discrete Values
2     16     382     38K     190K     1.7M     31M     2B

15.6s     43.1s     125.5s     226.6s     370.6s     690.5s     1329.4s     TO

10-Neuron Network Fair?

MiniACSIncome Details and  Verification Results.

§.§.§ Input Distribution

The Bayesian Network input distribution of MiniACSIncome has the network structure depicted in <ref>.
For using AGEP as a parent node, we summarised the age groups 17–34, 35–59, and 60–95.
This means that the conditional probability table of MAR does not have 78 entries for AGEP, but three, corresponding to the ranges 17–34, 35–59, and 60–95.

[minimum size=1.5cm, node distance=3cm,on grid]
[draw=black, circle] (RAC) at (0,0) RAC1P;
[draw=black, circle] (SEX) at (-3, 0) SEX;
[draw=black, circle] (AGEP) at (1.5, -2) AGEP;

[draw=black, circle] (SCHL) at (-1.5,-2)  SCHL;
[draw=black, circle] (MAR) at (0,-4) MAR;
[draw=black, circle] (RELP) at (1.5,-6) RELP;
[draw=black, circle] (COW) at (-1,-6) COW;
[draw=black, circle] (WKHP) at (-3.5,-6) WKHP;

[black,-Stealth] (SEX) – (SCHL);
[black,-Stealth] (RAC) – (SCHL);
[black,-Stealth] (SCHL) – (MAR);
[black,-Stealth] (AGEP) – (MAR);
[black,-Stealth] (MAR) – (RELP);
[black,-Stealth] (SCHL) – (COW);
[black,-Stealth] (MAR) – (COW);
[black,-Stealth] (SCHL) – (WKHP);
[black,-Stealth] (MAR) – (WKHP);

MiniACSIncome Bayesian Network Structure.

To fit the Bayesian Network, we walk the network from sources to sinks and fit each conditional distribution to match the empirical distribution of the data subset that matches the current condition.
For example, for fitting the conditional distribution of SCHL given SEX=1, RAC1P=1, we select the samples in the dataset having SEX=1 and RAC1P=1 and fit the conditional distribution of SCHL to match the empirical distribution of these samples.
We use categorical distributions for all categorical variables and WKHP.
For AGEP, we fit a mixture model of four truncated normal distributions that we discretise to integer values.

<Ref> depict the marginal distributions of 10 000 samples from the fitted Bayesian Network and the empirical marginal distributions of the MiniACSIncome-8 dataset.
<Ref> contains the correlation matrix of the same sample compared to the correlation matrix of MiniACSIncome-8.
The Bayesian Network approximates the empirical marginal distribution and the correlation structure of MiniACSIncome-8 reasonably well.

Note that for real-world fairness verification, the input distribution should not be fitted to the same dataset on which the neural network to verify is trained on.
Instead, the input distribution needs to be carefully constructed by domain experts.
For fairness audits, the input distribution could also be designed adversarially by a fairness auditing entity.

MiniACSIncome Bayesian Network — Marginal Distributions.
The generated data is sampled from the Bayesian Network, while the original data is the MiniACSIncome-8 dataset.

MiniACSIncome Bayesian Network — Correlation Matrix.
Population Model denotes the fitted Bayesian Network, while Training Data stands for the full MiniACSIncome-8 dataset.

§.§.§ Training

All MiniACSIncome neural networks are trained on a 56%/14%/30% split of the MiniACSIncome-$i$ dataset into training, validation, and testing data.
This split is identical for all MiniACSIncome-$i$ datasets.
All networks are trained using: the Adam optimiser <cit.>, cross entropy as loss function, no \(L_2\) regularisation (weight decay), \(\beta_1 = 0.9\), \(\beta_2 = 0.999\), \(\varepsilon = {10}^{-8}\) as suggested by KingmaBa2015, and a learning rate decay by $0.1$ after 2000 and 4000 iterations.
The learning rate and the number of training epochs are contained in <ref>.
For each network, we perform five random restarts and select the network with the lowest cross entropy on the validation data.

<Ref> contains the accuracy, precision, and recall for the overall test set, the persons with female sex in the test set, and the persons with male sex in the test set.
Additionally, the table contains whether a network satisfies the demographic parity fairness notion according to .

We used Optuna <cit.> for an initial exploration of the hyperparameter space but did not apply automatic hyperparameter optimisation to obtain the final training hyperparameters.

Dataset     Architecture     Learning Rate     # Epochs
MiniACSIncome-1     1$\times$10     0.0001     1

MiniACSIncome-2     1$\times$10     0.001      2

MiniACSIncome-3     1$\times$10     0.001      3

MiniACSIncome-4     1$\times$10     0.001      4

MiniACSIncome-5     1$\times$10     0.001      5

MiniACSIncome-6     1$\times$10     0.001      3

MiniACSIncome-7     1$\times$10     0.001      3

MiniACSIncome-8     1$\times$10     0.001      3

MiniACSIncome-4     1$\times$1000     0.001      2

MiniACSIncome-4     1$\times$2000     0.001      2

MiniACSIncome-4     1$\times$3000     0.001      2

MiniACSIncome-4     1$\times$4000     0.001      2

MiniACSIncome-4     1$\times$5000     0.001      2

MiniACSIncome-4     1$\times$6000     0.001      2

MiniACSIncome-4     1$\times$7000     0.001      2

MiniACSIncome-4     1$\times$8000     0.001      2

MiniACSIncome-4     1$\times$9000     0.001      2

MiniACSIncome-4     1$\times$10 000     0.001      2

MiniACSIncome-4     2$\times$10     0.001      2

MiniACSIncome-4     3$\times$10     0.001      2

MiniACSIncome-4     4$\times$10     0.001      2

MiniACSIncome-4     5$\times$10     0.001      2

MiniACSIncome-4     6$\times$10     0.001      2

MiniACSIncome-4     7$\times$10     0.001      2

MiniACSIncome-4     8$\times$10     0.001      2

MiniACSIncome-4     9$\times$10     0.001      2

MiniACSIncome-4     10$\times$10     0.001      2

MiniACSIncome Training Hyperparameters.

3cOverall     3cFemale     3cMale
Dataset     Net     A     P     R     A     P     R     A     P     R     Fair?
mACSI-1     1$\times$10           $ 57\% $     $ 44\% $     $ 63\% $     $ 72\% $         –       $  0\% $     $ 44\% $     $ 44\% $     $100\% $     \(\exfailure\)

mACSI-2     1$\times$10           $ 65\% $     $ 55\% $     $ 14\% $     $ 71\% $         –       $  0\% $     $ 58\% $     $ 55\% $     $ 23\% $     \(\exfailure\)

mACSI-3     1$\times$10           $ 73\% $     $ 67\% $     $ 48\% $     $ 76\% $     $ 65\% $     $ 35\% $     $ 69\% $     $ 68\% $     $ 55\% $     \(\exfailure\)

mACSI-4     1$\times$10           $ 75\% $     $ 68\% $     $ 60\% $     $ 78\% $     $ 66\% $     $ 48\% $     $ 72\% $     $ 69\% $     $ 66\% $     \(\exfailure\)

mACSI-5     1$\times$10           $ 76\% $     $ 69\% $     $ 63\% $     $ 79\% $     $ 65\% $     $ 53\% $     $ 74\% $     $ 71\% $     $ 69\% $     \(\exfailure\)

mACSI-6     1$\times$10           $ 77\% $     $ 69\% $     $ 63\% $     $ 79\% $     $ 65\% $     $ 55\% $     $ 74\% $     $ 72\% $     $ 68\% $     \(\exfailure\)

mACSI-7     1$\times$10           $ 77\% $     $ 71\% $     $ 64\% $     $ 79\% $     $ 67\% $     $ 50\% $     $ 76\% $     $ 72\% $     $ 72\% $     \(\exfailure\)

mACSI-8     1$\times$10           $ 78\% $     $ 70\% $     $ 68\% $     $ 80\% $     $ 67\% $     $ 56\% $     $ 76\% $     $ 71\% $     $ 75\% $     \(\exunknown\)

mACSI-4     1$\times$1000         $ 75\% $     $ 71\% $     $ 54\% $     $ 79\% $     $ 70\% $     $ 44\% $     $ 72\% $     $ 71\% $     $ 61\% $     \(\exfailure\)

mACSI-4     1$\times$2000         $ 75\% $     $ 65\% $     $ 67\% $     $ 78\% $     $ 61\% $     $ 60\% $     $ 72\% $     $ 67\% $     $ 71\% $     \(\exsuccess\)

mACSI-4     1$\times$3000         $ 74\% $     $ 63\% $     $ 71\% $     $ 77\% $     $ 58\% $     $ 67\% $     $ 72\% $     $ 66\% $     $ 73\% $     \(\exsuccess\)

mACSI-4     1$\times$4000         $ 75\% $     $ 71\% $     $ 55\% $     $ 79\% $     $ 69\% $     $ 44\% $     $ 72\% $     $ 71\% $     $ 61\% $     \(\exfailure\)

mACSI-4     1$\times$5000         $ 75\% $     $ 69\% $     $ 58\% $     $ 78\% $     $ 71\% $     $ 39\% $     $ 73\% $     $ 68\% $     $ 69\% $     \(\exfailure\)

mACSI-4     1$\times$6000         $ 75\% $     $ 69\% $     $ 59\% $     $ 79\% $     $ 69\% $     $ 44\% $     $ 72\% $     $ 68\% $     $ 67\% $     \(\exfailure\)

mACSI-4     1$\times$7000         $ 75\% $     $ 66\% $     $ 64\% $     $ 78\% $     $ 62\% $     $ 55\% $     $ 72\% $     $ 68\% $     $ 68\% $     \(\exfailure\)

mACSI-4     1$\times$8000         $ 75\% $     $ 66\% $     $ 64\% $     $ 78\% $     $ 63\% $     $ 55\% $     $ 72\% $     $ 68\% $     $ 69\% $     \(\exfailure\)

mACSI-4     1$\times$9000         $ 75\% $     $ 70\% $     $ 53\% $     $ 78\% $     $ 66\% $     $ 47\% $     $ 71\% $     $ 72\% $     $ 56\% $     \(\exfailure\)

mACSI-4     1$\times$10 000      $ 75\% $     $ 67\% $     $ 63\% $     $ 79\% $     $ 69\% $     $ 45\% $     $ 72\% $     $ 66\% $     $ 73\% $     \(\exfailure\)

mACSI-4     2$\times$10           $ 75\% $     $ 67\% $     $ 59\% $     $ 78\% $     $ 65\% $     $ 48\% $     $ 72\% $     $ 69\% $     $ 65\% $     \(\exfailure\)

mACSI-4     3$\times$10           $ 75\% $     $ 68\% $     $ 58\% $     $ 79\% $     $ 67\% $     $ 47\% $     $ 72\% $     $ 69\% $     $ 65\% $     \(\exfailure\)

mACSI-4     4$\times$10           $ 75\% $     $ 68\% $     $ 58\% $     $ 78\% $     $ 66\% $     $ 46\% $     $ 72\% $     $ 69\% $     $ 66\% $     \(\exfailure\)

mACSI-4     5$\times$10           $ 75\% $     $ 68\% $     $ 59\% $     $ 78\% $     $ 67\% $     $ 47\% $     $ 72\% $     $ 69\% $     $ 66\% $     \(\exfailure\)

mACSI-4     6$\times$10           $ 75\% $     $ 67\% $     $ 61\% $     $ 78\% $     $ 67\% $     $ 46\% $     $ 72\% $     $ 67\% $     $ 70\% $     \(\exfailure\)

mACSI-4     7$\times$10           $ 75\% $     $ 65\% $     $ 66\% $     $ 78\% $     $ 63\% $     $ 52\% $     $ 72\% $     $ 65\% $     $ 74\% $     \(\exfailure\)

mACSI-4     7$\times$10           $ 75\% $     $ 67\% $     $ 59\% $     $ 78\% $     $ 65\% $     $ 46\% $     $ 72\% $     $ 68\% $     $ 67\% $     \(\exfailure\)

mACSI-4     9$\times$10           $ 75\% $     $ 67\% $     $ 59\% $     $ 78\% $     $ 65\% $     $ 47\% $     $ 72\% $     $ 68\% $     $ 67\% $     \(\exfailure\)

mACSI-4     19$\times$10          $ 75\% $     $ 67\% $     $ 61\% $     $ 78\% $     $ 64\% $     $ 49\% $     $ 72\% $     $ 68\% $     $ 68\% $     \(\exfailure\)

MiniACSIncome Neural Networks.
The abbreviation mACSI-$i$ stands for MiniACSIncome-$i$.
We report the accuracy (A), precision (P), and recall (R) of each trained network (Net) for the whole test dataset (Overall), the persons with SEX=2 in the test set (Female), and the persons with SEX=1 in the test set (Male).
Additionally, we report whether the network satisfies the demographic parity fairness notion according to  (Fair?).

§ HEURISTICS FOR

This section contains motivation and additional implementation details for our   heuristic and compares the different branch and split selection heuristics that we introduce in <ref>.

§.§ BaBSB

Our  split selection heuristic is a variation of the  heuristic for non-probabilistic neural network verification of BunelLuTurkaslanEtAl2020.
One difference is that BunelLuTurkaslanEtAl2020 use the method of WongKolter2018 for estimating the improvement in bounds, while we use interval arithmetic.
Another difference is that while BunelLuTurkaslanEtAl2020 are mainly interested in lower bounds, we are equally interested in lower and upper bounds.
Therefore, we select the dimension \(d\) that yields the largest lower bound or smallest upper bound in any of the new branches, while BunelLuTurkaslanEtAl2020 select the dimension \(d\) with the largest lower bound among the smaller lower bound for the two branches originating from splitting \(d\).
We found this variant to be the most successful for our application.
BunelLuTurkaslanEtAl2020 discuss further variants.

We round all bounds to four decimal places to mitigate floating point issues.
If several dimensions yield equal improvements in bounds, we randomly select one of these dimensions.
Without this random tie-breaking, we might split a single dimension repeatedly if the interval arithmetic bounds are very loose.
We use a separate pseudo-random number generator with a fixed seed for this tie-breaking so that  remains entirely deterministic.

§.§ Experiments

We compare the branch selection heuristics  and  when running  on a selection of the benchmarks from <ref>.
In both cases, we use * to compute bounds on \(\pfun\) and use our  split selection heuristic.
We do not provide a detailed comparison of  and  since  does not allow  to solve even the least demanding FairSquare benchmarks within the time budget of 15 minutes.

Unlike the experiments in <ref>, we use recent hardware for comparing the different  heuristics.
Concretely, we use a compute server running Ubuntu 22.04 with an AMD Ryzen Threadripper 3960X CPU with 24 cores and 252 GB of memory (HW2).

We compare  and  on the FairSquare benchmarks, the ACAS Xu robustness benchmark, and MiniACSIncome-1 to MiniACSIncome-6.
We note that this set contains many less demanding instances.
<Ref> visualises the results.
has a slight advantage over , meaning that  can solve a few more instances and can solve a few other instances faster when using  compared to .

title=Solved Instances vs. Runtime,
ymin=0, ymax=1200,
xmin=0, xmax=170,
xlabel=# Solved Instances, ylabel=Runtime (s),

legend pos=south east,
legend cell align=left,
[Colors-A,very thick] table[x=Nr,y=Runtime,col sep=comma] data/Runtimes-prob.csv;

[Colors-B,very thick] table[x=Nr,y=Runtime,col sep=comma] data/Runtimes-prob-log-bounds.csv;

vs. . Lower is better.