We present a data-driven approach to the quantitative verification of probabilistic programs and stochastic dynamical models. Our approach leverages neural networks to compute tight and sound bounds for the probability that a stochastic process hits a target condition within finite time. This problem subsumes a variety of quantitative verification questions,
from the reachability and safety analysis of discrete-time stochastic dynamical models,
to the study of assertion-violation and termination analysis of probabilistic programs.
We rely on neural networks to represent supermartingale certificates that yield such probability bounds,
which we compute using a counterexample-guided inductive synthesis
we train the neural certificate while tightening the probability bound over samples of the state space using stochastic optimisation, and then we formally check the certificate's validity over every possible state using satisfiability modulo theories; if we receive a counterexample, we add it to our set of samples and repeat the loop until validity is confirmed.
We demonstrate on a diverse set of benchmarks
that, thanks to the expressive power of neural networks,
our method yields smaller or comparable probability bounds than existing symbolic methods in all cases,
and that our approach succeeds on models that are entirely beyond the reach of such alternative techniques.

§ INTRODUCTION

Probabilistic programs extend imperative programs with the ability to sample from probability distributions [22, 29, 33, 39],
which provides an expressive language to describe randomized algorithms, cryptographic protocols,
and Bayesian inference schemes.
Discrete-time stochastic dynamical models,
characterised by stochastic difference equations,
are a natural framework to describe auto-regressive time series, as well as sequential decision and planning problems in unknown environments.
A fundamental quantitative verification problem for probabilistic programs and stochastic dynamical models is the quantitative reachability question,
which amounts to finding the probability with which the system reaches a given target condition within a finite number of steps.
Reachability is at the core of a variety of other important quantitative verification questions as,
by selecting appropriate target conditions on the state space, we can express the probability that a probabilistic program
terminates or that it violates an assertion, as well as the probability that a stochastic dynamical model satisfies an invariant or remains within a set of safe configurations.

Quantitative reachability verification has been studied extensively using theories and algorithms built upon symbolic reasoning techniques,
such as quantitative calculi [38, 41],
probabilistic model checking algorithms [35, 24],
discrete abstractions of stochastic dynamical models [4, 49],
and the synthesis of supermartingale-like certificates [15, 17, 21, 18].
Among the latter class, a method to provide a sound upper bound for the reachability probability of a system
is to synthesise a supermartingale function that maps every reachable state to a non-negative real,
whose value is never smaller than 1 inside the target condition,
and such that it never increases in expectation as the system evolves outside of the target.
This is referred to as a non-negative repulsing supermartingale or stochastic invariant indicator in the literature [48, 19, 21],
and its output over a given state provides an upper bound for the probability that the system reaches the target condition from that state.
Symbolic methods for the synthesis of such certificates assume that the supermartingale function,
as well as its post-expectation, which depends on the model constraints and distributions,
are both in linear or polynomial form. This poses syntactic restrictions to their applicability.

Data-driven and counterexample-guided inductive synthesis (CEGIS) procedures,
combined with machine learning approaches that leverage neural networks to represent certificates,
have shown great promise in mitigating the aforementioned limitation [16, 1, 28, 43, 37, 3, 36].
In particular, neural-based CEGIS decouples the task of guessing a certificate from that of checking its validity,
delegating the guessing task to efficient machine learning algorithms that leverage the expressive power of neural networks,
while confining symbolic reasoning to the checking part of the task,
which is computationally easier to solve in isolation than the entire synthesis problem.
CEGIS has been applied to the synthesis of neural supermartingales for the almost-sure termination of
probabilistic programs [3],
as well as their counterpart for stochastic dynamical models with applications to qualitative queries such as almost-sure
safety and stability [37, 36].

In this paper, we present theory, methods, and an extensive experimental evaluation to demonstrate the efficacy and flexibility of neural supermartingales to solve quantitative verification questions for probabilistic programs and stochastic dynamical models, using machine learning combined with satisfiability modulo theories (SMT) technologies.

Theory We adapt the theory of non-negative repulsing supermartingales to leverage neural networks
as representations of supermartingale functions for quantitative verification. Unlike previous deductive methods
which need deterministic invariants and impose restrictions on the models under study,
our version entirely relies on neural architectures and SMT solving to guarantee soundness,
without requiring such assumptions.

Methods We present a CEGIS-based approach to train neural supermartingale functions that minimise an upper bound for the reachability  probability over sample points from the state space, and check their validity over every possible state using SMT solving.
We present a program-agnostic approach that relies on state samples in the training phase, and also a novel
program-aware approach that embeds model information in the loss function to enhance the effectiveness of stochastic optimisation.

Experiments We build a prototype implementation and compare the efficacy of our method with the
state of the art in synthesis of linear supermartingales using symbolic reasoning [48].
We show that our program-aware approach computes tighter or comparable
probability bounds than symbolic reasoning on existing benchmarks, while our program-agnostic approach matches it in over half of the instances.
We additionally demonstrate that both our approaches can handle models beyond reach of purely symbolic methods.

§ PROBABILISTIC PROGRAMS AND STOCHASTIC DYNAMICAL MODELS

\begin{align*}
v &\in \mbox{Vars} & \text{(variables)}\\
N &\in \Real & \text{(numerals)}
\\E &::= v \mid N \mid \texttt{-}E \mid E+E \mid E-E \mid E*E \mid \dots
&\text{(arithmetic expressions)}
\\P &::= {\tt Bernoulli}\texttt{(}\,E\,\texttt{)} \mid
{\tt Gaussian}\texttt{(}\,E\texttt{, }E\,\texttt{)} \mid \dots
&\text{(probability distributions)}
\\B &::= {\tt true} \mid {\tt !}B \mid B \,\texttt{\&\&}\, B \mid B\texttt{||}B  \mid E\,\texttt{==}\,E \mid E\,\texttt{<}\,E \mid \dots & \text{(Boolean expressions)}
\\C &::= \keyword{skip} &\text{(update commands)}
\\&\quad\mid v\,\texttt{=}\,E&\text{(deterministic assignment)}
\\&\quad\mid v\,\probeq\,P&\text{(probabilistic assignment)}
\\&\quad\mid C~{\tt ;}~C &\text{(sequential composition)}
\\&\quad\mid \keyword{if}~B~\keyword{then}~C~\keyword{else}~C~\keyword{fi}&\text{(conditional composition)}
\end{align*}

Grammar for update commands, Boolean and arithmetic expressions.

Probabilistic programs are computer programs whose execution is determined by random variables,
and stochastic dynamical models describe discrete-time dynamical systems with probabilistic behaviour.
The earlier enjoy the flexibility of imperative programming constructs and are used to describe randomised algorithms,
and the latter are expressed as stochastic difference equations and are used to describe probabilistic systems that evolve
over infinite time. The semantics of both can be described in terms of  stochastic processes and,
for this reason, verification questions for both  can be solved with similar techniques.

The syntax of our modeling framework uses imperative constructs from probabilistic programs and defines
executions over infinite time as dynamical systems.
Specifically, we consider programs that operate over an ordered set of $n$ real-valued variables, denoted by $\Vars$, and update their values through the repeated execution of a command $C$ whose grammar is described in Figure <ref>.
Under this definition, a state of the
system is an $n$-dimensional vector $s \in \Real^n$ that assigns a value to each variable symbol. The update command $C$
defines an update function $f \colon \Real^n \times [0,1]^m \to \Real^n$, where $m$ is the number of syntactic probabilistic assignment statements occurring in $C$, which maps the current state and $m$ random variables uniformly distributed in $[0, 1]$ into the next state.
Conceptually, within $f$, each random variable is mapped into its respective distribution by applying the appropriate inverse transformation.
Altogether, our probabilistic program defines a stochastic process, whose
behavior is determined by the following stochastic difference equation:

\begin{equation}
s_{t+1} = f(s_t, r_t), \quad r_t \sim \mathbb{U}^m,
\label{eq:sys}
\end{equation}

where $r_t$ is an $m$-dimensional random input sampled at time $t$
from the uniform distribution $\mathbb{U}^m$ over the $m$-dimensional hypercube $[0,1]^m$.
The initial state $s_0$ is either given as a deterministic assignment
to constant values, or is non-deterministically chosen from a set of initial
conditions $S_0$ characterized by a Boolean expression.
This setting can be seen as an assertion about the initial conditions followed by the probabilistic program
As we show in <ref>, this allows us to characterise verification
questions such as termination, non-termination, invariance and
assertion-violation for while loops with general guards,
as well as reachability and safety verification questions for dynamical models with general stochastic disturbances.

The semantics of our model is defined as a stochastic process induced by the Markov chain over the probability space of infinite words of random samples.
This is defined by the probability space triple $(\Omega, \mathcal{F}, \mathbb{P})$ [14, 30], where

* $\Omega$ is the set of infinite sequences $\left([0, 1]^m\right)^\omega$ of $m$-dimensional tuples of values in $[0,1]$,

* $\mathcal{F}$ is the extension of the Borel $\sigma$-algebra over the unit interval $\mathcal B([0, 1])$ to $\Omega$,

* $\mathbb{P}$ is the extension of the Lebesgue measure on $[0, 1]$ to $\Omega$.

Every initialisation of the system on state $s \in \Real^n$ induces a stochastic
process $\{ X_t^{s} \}_{t \in \Nat}$ over the state space $\Real^n$.
Let $\omega = r_0r_1r_2\dots$ be an infinite sequence of random samples in $[0,1]^m$,
then the stochastic process is defined by the sequence of random variables

\begin{equation}\label{eqn:stochprocess}
X_{t + 1}^{s}(\omega) =
f(X_{t}^{s}(\omega), r_t),
\qquad X_0^{s}(\omega) = s.
\end{equation}

This defines the natural filtration $\{ \mathcal{F}_t \}_{t \in \Nat}$,
which is the smallest filtration to which the stochastic process $X^s_t$ is adapted.
In other words [32], this can be seen as another Markov chain with state space $\Real^n$ and transition kernel

\begin{equation}
T(s, S') = \Leb\left( \left\{r \in [0, 1]^m \mid f(s,r) \in S'\right\} \right),
\end{equation}

where $S'$ is a Borel measurable subset of $\Real^n$
and Leb refers to the Lebesgue measure of a measurable subset of $[0, 1]^m$. In other words, kernel $T$ denotes the probability to transition from state $s$ into a set of states $S'$.  The transition kernel also defines
the post-expectation operator $\mathbb X$, also known as the next-time operator <cit.>.
$\mathbb X$ can be applied to an arbitrary Borel-measurable function $h \colon \Real^n \to \Real$
defining the post-expectation of $h$, denoted by $\mathbb{X}[h]$ and defined as the following function over states:

\begin{equation}\label{eqn:postex}
\mathbb{X}[h](s)= \int %_{\Real^n}
h(s')~T(s, \mathrm{d}s').
\end{equation}

This represents the expected value of $h$ evaluated at the next state, given the current state being $s$.
Computing the symbolic representation of a post-expectation for probabilistic programs and stochastic dynamical models
is a core problem in probabilistic verification [26, 27]. Indeed, our theoretical framework builds upon the post-expectation (cf. Eq. (<ref>)),
and our verification procedure uses a symbolic representation of the post-expectation (cf. <ref>) as well as our program-aware learning procedure, while our program-agnostic learning procedure approximates it statistically (cf. <ref>).

We remark that our model encompasses general probabilistic program loops as well as stochastic dynamical systems
with general disturbances. For example, a probabilistic loop with guard condition $B$ and body $C$ in the form

\begin{align}
\keyword{while}~B~\keyword{do}~C~\keyword{od} \label{eq:loop}
\end{align}

can be expressed as a loop with guard $\keyword{true}$ and
body $\keyword{if}~B~\keyword{then}~C~\keyword{else}~\keyword{skip}~\keyword{fi}$.
This expresses the fact that, after termination, the program will stay on the terminal state indefinitely.
Also, discrete-time stochastic difference equations with a nonlinear vector field $g$ and time-invariant input disturbance with an arbitrary distribution ${\cal W}$ in the form

\begin{equation}
s_{t+1} = g(s_t, w_t), \quad w_t \sim \mathcal{W}. \label{eq:dynsys}
\end{equation}

comply with our model. It is sufficient to derive $w_t$ with an appropriate
inverse transformation from the uniform distribution and embed it in $f$.
Our model is even more general, as it
encompasses state-dependent distributions,
whose parameters depend on the state and may depend on other distributions and thus define joint, multi-variate
and hierarchically-structured distributions. Notably, our model comprises both continuous and discrete probability distributions and is able to model discrete-time stochastic hybrid systems.

§ QUANTITATIVE REACHABILITY VERIFICATION OF PROBABILISTIC MODELS

Quantitative verification treats the question of providing a quantity
for which a system satisfies a property as opposed to providing a definite positive or negative answer. In fact, it is sometimes too conservative or inappropriate to demand a definite outcome to
a formal verification question. For instance, a system whose behaviour is probabilistic may violate a specification on rare corner cases and yet satisfy it with a probability that is deemed acceptable for the application domain. We address the quantitative verification of probabilistic systems, which is the problem of computing the probability for which a system satisfies a specification [34, 50, 10, 9].

We consider the quantitative reachability verification question, which as we show below, is at the core of a variety of quantitative verification questions for probabilistic programs and stochastic dynamical models. Henceforth, we use $\mathbf{1}_S$
to denote the indicator function of set $S$, i.e.,
$\mathbf{1}_S(s) = 1$ if $s \in S$ and $\mathbf{1}_S(s) = 0$ if $s \not\in S$; we also use $\lambda x.M$ to denote the anonymous function that takes an argument $x$ and evaluates the expression $M$ to produce its result.
We now characterise the probability that a stochastic process reaches a Borel measurable target set $A$ in exactly time $t$, in at most time $t$, and in any finite time.

Let the event that a stochastic process $\{ X^s_t\}_{t \in \Nat}$
over state space $\Real^n$
initialised in state $s \in \Real^n$
reaches a target set $A \in \mathcal B(\Real^n)$
in exactly time $t \in \Nat$ be

\begin{equation}\label{eqn:reach1}
\Reach_t^{s}(A) = \left\{ \omega \in \Omega \mid X_0^{s}(\omega) \not\in A , \ldots, X_{t-1}^{s}(\omega) \not\in A , X_{t}^{s}(\omega) \in A \right\},
\end{equation}

with $\Reach_0^{s}(A) = \Omega$ if $s \in A$, and $\Reach_0^{s}(A) = \emptyset$ if $s \notin A$.
Then, $\Reach_t^{s}(A)$ is measurable and its probability measure can be expressed as follows:

\begin{align}
\mathbb{P}[\Reach_{{t+1}}^{s}(A)]
\mathbf{1}_{\mathbb{R}^n \setminus A}(s)
\cdot
\mathbb{X}[\lambda s' \mathpunct{.}
\mathbb{P}[\Reach_{{t}}^{s'}(A)]
\label{eqn:standard-result-reach}\\
\mathbb{P}[\Reach_{0}^{s}(A)]
&= \mathbf{1}_A(s).
\label{eqn:standard-result-zero}
\end{align}

Let the event that a stochastic process $\{ X^s_t \}_{t \in \Nat}$
over state space $\Real^n$
initialised in state $s \in \Real^n$
reaches a target set $A \in \mathcal B(\Real^n)$
in at most time $t \in \Nat$

\begin{equation}\label{eqn:reach-union}
\Reach_{\leq t}^{s}(A) = \cup \{~\Reach_{i}^{s}(A) \colon 0 \leq i \leq t~\},
\end{equation}

Then, $\Reach_{\leq t}^{s}(A)$ is measurable and its probability measure can be expressed as follows:

\begin{align}
\mathbb{P}[\Reach_{\leq t + 1}^{s}(A)]
\mathbf{1}_A(s)
\mathbf{1}_{\Real^n \setminus A}(s)
\cdot
\mathbb{X}[\lambda s' \mathpunct{.}
\mathbb{P}[\Reach_{\leq t}^{s'}(A)]](s),
\label{eqn:reachafterzero}\\
\mathbb{P}[
\Reach_{\leq 0}^{s}(A)
] &= \mathbf{1}_A(s).
\label{eqn:reachbeforezero}
\end{align}

Let the event that a stochastic process $\{ X^s_t \}_{t \in \Nat}$
over state space $\Real^n$
initialised in state $s \in \Real^n$
reaches a target set $A \in \mathcal B(\Real^n)$
in finite time be

\begin{equation}\label{eqn:reachfin}
\Reach_{\mathrm{fin}}^s(A) = \cup \left\{\Reach_{ t}^{s}(A) \colon t \in \Nat\right\} =
\cup \left\{\Reach_{\leq t}^{s}(A) \colon t \in \Nat\right\}.
\end{equation}

Then, $\Reach_{\mathrm{fin}}^s(A)$ is measurable and its probability measure (which we refer to as the reachability probability of the target set) can be expressed as follows:

\begin{equation}\label{eqn:reachfinlim}
\mathbb P [\Reach_{\mathrm{fin}}^s(A)] = \lim_{t \to \infty}\mathbb P [\Reach_{\leq t}^s(A)].
\end{equation}

The lemmata above follow from measure-theoretic results in probabilistic verification (proofs are provided in <ref>), and underpin the formal characterisation of quantitative probabilistic reachability in the next section.
Our method leverages neural networks to compute an upper bound
for the reachability probability with respect to a given target set (cf. <ref>).
Notice that an appropriate choice of target set allows us to express a variety of other quantitative verification questions,
for which our method can provide an upper or lower bound.
Below, by providing a suitable choice for the target set, we show how important quantitative verification questions for probabilistic programs and stochastic dynamical models can be characterised as instances of quantitative reachability.

Termination Analysis
Let $G \in \mathcal B (\Real^n)$ be the guard set of a probabilistic while loop as in Eq. (<ref>), i.e., the set of states for which
that guard condition evaluates to true.
The event that the loop terminates from initial state $s_0$ is $\Reach_{\mathrm{fin}}^{s_0}(\Real^n \setminus G)$.
Our method computes an upper bound for the probability that the loop terminates or, dually, it computes a lower bound for the probability of non-termination. Notably, when this lower bound is greater than 0, then almost-sure termination is refuted.

Assertion-violation Analysis
Let $G \in \mathcal B (\Real^n)$ be the guard set of a
probabilistic while loop and $A \in \mathcal B (\Real^n)$ be the satisfying set of an assertion placed at the beginning of the
loop body. Given initial state $s_0$, the event that the assertion is eventually violated
is $\Reach_{\mathrm{fin}}^{s_0}(G \setminus A)$.
Our method computes an upper bound for the probability of assertion violation or, dually, a lower bound for its satisfaction. Note that  assertions in other positions of the body can be modelled similarly by using additional Boolean variables.

Safety Verification
Let $B \in \mathcal B (\Real^n)$ be a set of undesirable states in a stochastic dynamical model. The event that the system is safe when initialised in $s_0$, i.e., it never reaches an undesirable state,
is given by
$\Omega \setminus \Reach_{\mathrm{fin}}^{s_0}(B)$.
Our method computes a lower bound for the probability that the system is safe, which is $1 - \mathbb P [\Reach_{\mathrm{fin}}^{s_0}(B)]$.

Invariant Verification Let $I \in \mathcal B (\Real^n)$ be a candidate invariant set. The event that $I$ is invariant when the system is initialised in $s_0$ is $\Omega \setminus \Reach_{\mathrm{fin}}^{s_0}(\Real^n \setminus I)$.
Our method computes a lower bound for the probability
that set $I$ is invariant, which is $1 - \mathbb P [\Reach_{\mathrm{fin}}^{s_0}(\Real^n \setminus I)] $. Note that if $s_0 \notin I$, this definition yields a trivial lower bound of zero for the probability of invariance.

§ NEURAL SUPERMARTINGALES FOR QUANTITATIVE VERIFICATION

Supermartingale certificates provide a flexible and powerful
theoretical framework for the formal verification of
probabilistic models with infinite state spaces.
Not only have supermartingales been applied to qualitative questions,
such as almost-sure termination analysis of probabilistic programs
and almost-sure stability analysis of stochastic dynamical models,
but also to quantitative reachability verification, as in this work.
Specifically, this is enabled by the theory of non-negative repulsing supermartingales and of stochastic invariants [48, 19, 21].

Our method builds upon the theory of non-negative repulsing supermartingales and stochastic invariant indicators
and adapts it to take advantage of the expressive power of neural networks and the
flexibility of machine learning and counterexample-guided inductive synthesis algorithms. Our method hinges on the following theorem, whose proof we provide in <ref>.

Let $\mathbb{X}$ be the post-expectation operator of a stochastic process
over state space $\Real^n$ and let $A \in \mathcal{B}(\Real^n)$ be a target set. Let $V \colon \Real^n \to \Real_{\geq 0} \cup \{\infty\}$ be a non-negative function that satisfies the following two conditions:

(indicating condition)       ∀s ∈A V(s) ≥1,

(non-increasing condition)       ∀s ∉A 𝕏[V](s) ≤V(s),      aaaaaaaa

Then, for every state $s \in \Real^n$, it holds that $V(s) \geq \mathbb P [\Reach_{\mathrm{fin}}^{s}(A)]$.

As we show in detail in <ref>, we compute a supermartingale that satisfies the two criteria (<ref>) and (<ref>),
while also minimising its output over the initial state $s_0$. When the initial state is chosen nondeterministically from a set $S_0$, we instead minimise the output over all states in the set. As a consequence, the maximum of $V$ over $S_0$ is a sound upper bound for the reachability probability.

In this work, we template $V$ as a neural network with $n$ input neurons, $l$ hidden layers with respectively
$h_1, \dots h_l$ neurons in each hidden layer, and 1 output neuron.
We guarantee a-priori that the function's output will be non-negative over the entire domain using the following architecture:

\begin{equation}
\label{eq:net}
V(x) = \operatorname{sum} \left( \sigma_l \circ \dots \circ \sigma_1(x) \right), \quad
\sigma_i(z) =
(\relu(\mathbf{w}^T_{i,1} z + \mathbf{b}_{i,1}), \dots, \relu(\mathbf{w}^T_{i,h_i} z + \mathbf{b}_{i,h_i}))
\end{equation}

where $\operatorname{sum}(z_{l, 1}, \ldots, z_{l, h_l}) = \sum_{k = 1}^{h_l} z_{l, k}$ and $\mathbf{w}^T_{i,j} \in \mathbb{Q}^{h_{i-1}}$ (defining $h_0 = n$, the number of input neurons) and
$\mathbf{b}_{i,j} \in \mathbb{Q}$ are respectively the weight vector and the bias parameter for the inputs to neuron $j$ at layer $i$.
Then, our method trains the neural network to satisfy the two criteria (<ref>) and (<ref>),
while also minimising its output on the initial condition.

while x > 0 do
assert(x <= 10);
p   Bernoulli(0.5);
if p == 1 then
x -= 2
x += 1

[outer sep=0, inner sep=0]
[axis x line=center,
axis y line=center,
xtick=-2, -1, 0,1,...,11,
x label style=anchor=north, below=5mm, left=17mm,
y label style=rotate=90,anchor=south, left=20mm, above=1mm,
xlabel=Initial value of x,

[black!60, very thick] table ;
[orange, very thick, domain=-3:11] x/12+1/6;
[blue,very thick] coordinates (-3, 0.00) (-2,0.00) (0,0.04) (1,0.06) (2,0.09) (3,0.12) (4,0.16) (5,0.21) (6,0.27) (7,0.36) (8,0.5) (9,0.7) (10,1) (10.5, 1.15);
[very thick, magenta] coordinates (-3, 0.01) (-2, 0.01) (0, 0.069) (5, 0.39) (10, 1) (10.5, 1.061);

at (2.9,2.2) $V^{\text{lin}}$;
at (3.0,1.55) $V^{\text{neur}}_{3}$;
at (3.0,1) $V^{\text{neur}}_{12}$;
at (4.5, 0.8)  $P_\text{reach}$ ;

Comparison between linear and neural supermartingale functions in tightness of bounds for the assertion violation probability of the program repulse, shown on the left. On the right, the function ${V^\text{lin}}$ indicates the tightest linear supermartingale (under the restriction that x is greater than $-2$). The functions ${V^\text{neur}_3}$ and ${V^\text{neur}_{12}}$ indicate single-layer neural supermartingales with 3 and 12 neurons respectively. The piecewise constant function $P_\text{reach}$ is the true probability of assertion violation, and indicates the ideal lower bound for the value of any supermartingale function.

Using neural networks as templates of supermartingale functions
introduces non-trivial advantages with respect to symbolic methods
for the synthesis of supermartingales.
Firstly, neural supermartingales are able to better approximate the true probability of reachability and thus attain tighter upper bounds on it.
Secondly, symbolic methods require deterministic invariants that overapproximate the set of reachable states, and suitably restricts the domain of the template. Figure <ref> illustrates using an example the advantages of neural certificates with respect to linear supermartingales synthesised using Farkas' lemma.
This example shows that increasing the number of neurons provides greater flexibility and allows the certificate to more tightly approximate the true probability.
Moreover, this example shows that linear supermartingales require their domain to be restricted with an appropriate deterministic invariant.
Notably, symbolic methods for the synthesis of polynomial supermartingales based on Putinar's Positivstellensatz also require compact deterministic invariants to be provided [19]. By contrast, neural supermartingales (whose output is always non-negative) achieve the same result while
relaxing the requirement of providing an invariant beforehand.

§ DATA-DRIVEN SYNTHESIS OF NEURAL SUPERMARTINGALES

Our approach to synthesising neural supermartingales for quantitative verification utilises a counterexample-guided inductive synthesis (CEGIS) procedure [46, 47] (cf. <ref>). This procedure consists of two components,
a learner and a verifier, that work in opposition to each other.
On the one hand, the learner seeks to synthesise a candidate supermartingale that meets the desired specification (cf. Eq. (<ref>)) over a finite set of samples from the state space, while simultaneously optimising the tightness of the probability bound. On the other hand, the verifier seeks to disprove the validity of this candidate by searching for counterexamples, i.e., instances where the desired specification is invalidated, over the entire state space.
If the verifier shows that no such counterexample exists, then the desired specification is met by the supermartingale and the procedure provides a sound probability bound for the reachability probability of interest, together with a neural supermartingale to certify it.

(input) at (-2, 2) ;
(learner) [draw, rectangle, rounded corners, minimum width=1.5cm, minimum height =0.75cm]
at (0,2) Learner;

(verif) [draw, rectangle, rounded corners, minimum width=1.5cm, minimum height =0.75cm]
at (3,2) Verifier;

(out) [label=[align=left]$V, {p}$] at (5,2)  ;

[-Stealth[scale=1.5]]  (input) edge node[below, xshift=-0.2cm] $f$, $S_0$, $A$ (learner);

[-Stealth[scale=1.5], bend left] (verif) edge node[below] $d_\text{cex}$ (learner);

[-Stealth[scale=1.5], bend left] (learner) edge node[above] Candidate $V$ (verif);

[-Stealth[scale=1.5]] (verif) edge (out);

Overview of the counterexample-guided inductive synthesis procedure used to synthesise neural supermartingales for quantitative reachability verification. Inputs to the procedure are a probabilistic program $f$, a set of initial states $S_0$, and a target set $A$. The procedure outputs a valid neural supermartingale $V$ and a probability bound $p$.

§.§ Training of Neural Supermartingales From Samples

Our neural supermartingale for quantitative reachability verification consists of a neural network with $\relu$ activation functions, with an arbitrary number of hidden layers (cf. <ref>). We train this neural network using gradient descent over a finite set $D = \{ d^{(1)}, \dots , d^{(m)}\} \subseteq \mathbb{R}^n$ of states $d^{(i)}$ sampled over the state space $\Real^n$.
Initially, we sample uniformly within a bounded hyper-rectangle of $\mathbb{R}^n$, a technique that scales effectively to high dimensional state spaces and efficiently populates the initial dataset of state samples. Then, we construct a loss function that
guides gradient descent to optimise the parameters (weight and biases)
of the neural network $V$ to satisfy the specification set out in Eq. (<ref>) while also minimising the probability bound.
We define a loss function $\mathcal{L}(D)$ that consists of three terms:

\begin{equation}
\label{eq:loss}
\mathcal{L}(D) = \beta_{1}\mathcal{L}_\indicating(D) + \beta_{2} \mathcal{L}_\nonincreasing(D) + \beta_3\mathcal{L}_{\bound}(D).
\end{equation}

Components $\mathcal{L}_\indicating$ and $\mathcal{L}_\nonincreasing$ are responsible for encouraging satisfaction of the conditions in Eq. (<ref>), while the component $\mathcal{L}_{\bound}$ is responsible for tightening the probability bound. The parameters of this optimisation problem are the parameters of the neural network, which are initialised randomly. The dataset consists of the state samples, initially sampled randomly, and generated from counterexamples in subsequent CEGIS iterations (cf. <ref>). The coefficients $\beta_{1}, \beta_2$ and $\beta_3$ denote scale factors for each term, which we choose according to the priority that we want to assign to each condition (cf. <ref>).

First, consider the condition in Eq. (<ref>), which we refer to as the indicating condition. For this, we use the following loss function:

(indicating loss)
ℒ_(D) = 𝔼_d ∈D ∩A[(1 - V(d))].

This adds a penalty for states $d \in D$ lying inside the target condition $A$, at which $V$ fails to satisfy the indicating condition, whilst ignoring any states where $V$ satisfies it. We average this per-state penalty across all states in $D \cap A$.

We next consider the non-increasing condition in Eq. (<ref>), for which we use the following loss term:

(non-increasing loss)
ℒ_(D) = 𝔼_d ∈D ∖A[(𝕏[V](d) - V(d))].

This penalises states $d \in D$ lying outside of the target set $A$, at which $V$ fails to satisfy the non-increasing condition.
Notably, this component is defined in terms of the post-expectation $\mathbb{X}[V]$ of our supermartingale. To embed this expression in
our loss function we consider two alternative approaches, which we call program-aware and program-agnostic.

Program-aware Approach The program-aware approach uses the source code of the program to generate a symbolic expression for the post-expectation of $V$. For this purpose, we exploit the symbolic inference algorithm introduced by the tool PSI [26, 27], along with a symbolic representation of $V$. We construct a probabilistic program which represents the evaluation of $V$ on the state resulting after the execution of the update function $f$. The expected value of this program is precisely $\mathbb{X}[V]$. This results in a symbolic expression that is a function of the program state and parameters of $V$. In the non-increasing loss, states are instantiated to elements of $D \setminus A$, while the parameters are left as free-variables that the gradient descent engine differentiates with respect to.

Program-agnostic Approach
The program-agnostic approach provides an alternative formulation of the non-increasing loss term that does not require symbolic reasoning. Instead of leveraging the program's source code, it only requires the ability to execute it. For this, we utilise a Monte Carlo scheme to estimate the post-expectation.
For each state $d$ in our dataset $D$, to obtain an estimate of $\mathbb{X}[V](d)$ we sample a number $m^\prime$ of successor states $D^\prime = \{ {d^\prime}^{(k)}~:~1 \leq k \leq m^\prime\}$. Each successor state ${d^\prime} ^{(k)}$ is sampled by executing the program's update function $f$ (cf. Eq. (<ref>)) at state $d$. Then $\mathbb{X}[V](d)$ is estimated as
\approx \mathbb E_{d^\prime \in D'}[V(d^\prime)]$
which is the average of $V$ over $D^\prime$.
Even though this is an approximation, we emphasise that this does not affect the soundness of our scheme, which is ensured by the verifier.

Finally, we introduce a tightness criterion that minimises the probability upper bound:

(minimisation loss)

ℒ_(D) = 𝔼_d ∈D ∩S_0[V(d)].

This term encourages $V$ to take smaller values over the set of initial states $S_0$. Recall that the smaller the probability upper bound, the closer it is to the true value.

The loss function is provided to the gradient descent optimiser, whose performance benefits from a smooth objective.
For this reason, to improve the performance of our learner we replace every ReLU with a smooth approximation, Softplus, which takes the form $\softplus(s) = \log(1 + \exp(s))$. Additionally, we improve the approximation of Softplus to ReLU at small values over the interval $[0, 1]$ by re-scaling $V$. This means that we modify the bounding condition to require that $V(x) \geq \alpha$ over the target set $A$, for some large $\alpha > 1$. In other words, we modify the indicating loss to $\mathcal{L}_\indicating(V) = \mathbb E_{d \in D \cap A}[\relu(\alpha - V(d))]$.
We remark that while $\softplus$ is used as the activation function in the learning stage, for the verification stage we instead employ $\relu$ activation functions, ensuring soundness of the generated neural supermartingale. We remark that this has no effect on the soundness of our approach, which is ultimately guaranteed by the verifier.

§.§ Verification of Neural Supermartingales Using SMT Solvers

The purpose of the verification stage is to check that the neural supermartingale meets the requirements of Eq. (<ref>) over the entire state space $\mathbb{R}^n$, and if that is determined to be the case to furthermore obtain a sound upper bound on the reachability probability.
We achieve this by constructing a suitable formula in first-order logic, Eq. (<ref>), and use SMT solving to decide its validity, or equivalently, to decide the satisfiability of its negation, Eq. (<ref>).

The conditions pertaining to the validity of the neural supermartingale, given in (<ref>) and (<ref>), are encoded by the formulas
and $\varphi_\nonincreasing$. We also note that for a constant $p \in [0, 1]$ to be a sound upper bound on the reachability probability, it is sufficient to require that $p$ is an upper bound on the neural supermartingale's value over the set of initial states $S_0$ (cf. <ref>), which is expressed by the formula $\varphi_\verifbound$. A suitable choice for the bound $p$ is determined by a binary search over the interval $[0, 1]$.

\begin{equation}
\label{eq:val-q}
\forall s \in \Real^n: \underbrace{(s \in A \rightarrow V(s) \geq 1 )}_{\varphi_\indicating} \wedge \underbrace{(s \notin A \rightarrow \mathbb{X}[V](s) \leq V(s))}_{\varphi_\nonincreasing}
\wedge \underbrace{(s \in {S}_0 \rightarrow V(s) < p)}_{\varphi_\verifbound}.
\end{equation}

Here, $V(s)$ is a symbolic encoding of the candidate neural supermartingale proposed by the learner (<ref>), and $S_0$ and $A$ are defined by Boolean predicates over program variables, all of which (in our setting of networks composed from ReLU activations) are expressible using expressions and constraints in non-linear real arithmetic. For this reason, we use Z3 as our SMT solver [23].

The SMT solver is provided with the negation of Eq. (<ref>), namely

\begin{equation}
\label{eq:smt-spec}
\exists s \in \Real^n: \underbrace{(s \in A \wedge V(s) < 1 )}_{\lnot \varphi_\indicating} \vee \underbrace{(s \notin A \wedge \mathbb{X}[V](s) > V(s))}_{\lnot \varphi_\nonincreasing}
\vee \underbrace{(s \in {S}_0 \wedge V(s) \geq p)}_{\lnot\varphi_\verifbound},
\end{equation}

and decides its satisfiability, seeking an assignment $d_\text{cex}$ of $s$ that is a counterexample to the neural supermartingale's validity, for which any of $\lnot\varphi_\indicating$, $\lnot\varphi_\nonincreasing$ and $\lnot\varphi_\verifbound$ are satisfied. If no counterexample is found, this certifies the validity of the neural supermartingale.
Alternatively, if a counterexample $d_\text{cex}$ is found, it is added to the data set $D$, for the synthesis to incrementally resume.

§ EXPERIMENTAL EVALUATION

Results comparing neural supermartingales with Farkas Lemma for different benchmarks. Here, $p$ is the average probability bound generated by the certificate; success ratio is the number of successful experiments, out of 10 repeats, generated by CEGIS with neural supermartingale; `-' means no result was obtained. We also denote the architecture of the network: $(h_1, h_2)$ denotes a network with 2 hidden layers consisting of $h_1$ and $h_2$ neurons respectively.

0.86Benchmark     Farkas'     4cQuantitative Neural Certificates     Network
Lemma     2cProgram-Agnostic     2cProgram-Aware
(lr)3-4 (lr)5-6
$p$     Success Ratio      $p$     Success Ratio

persist_2d            -                 $\leq 0.1026$     0.9     $\leq 0.1175$     0.9      (3, 1)

faulty_marbles        -                 $\leq 0.0739$     0.9     $\leq 0.0649$     0.8     3

faulty_unreliable     -                 $\leq 0.0553$     0.9     $\leq 0.0536$     1.0     3

faulty_regions        -                 $\leq 0.0473$     0.9     $\leq 0.0411$     0.9     (3, 1)
cliff_crossing        $\leq 0.4546$     $\leq 0.0553$     0.9     $\leq 0.0591$     0.8     4

repulse100                $\leq 0.0991$     $\leq 0.0288$     1.0     $\leq 0.0268$     1.0     3

repulse100_uniform       $\leq 0.0991$     $\leq 0.0344$     1.0     -                 -       2

repulse100_2d            $\leq 0.0991$     $\leq 0.0568$     1.0     $\leq 0.0541$     1.0     3

faulty_varying        $\leq 0.1819$     $\leq 0.0864$     1.0     $\leq 0.0865$     1.0     2

faulty_concave        $\leq 0.1819$     $\leq 0.1399$     1.0     $\leq 0.1356$     0.9     (3, 1)
fixed_loop            $\leq 0.0091$     $\leq 0.0095$     1.0     $\leq 0.0094$     1.0     1

faulty_loop           $\leq 0.0181$     $\leq 0.0195$     1.0     $\leq 0.0184$     1.0     1

faulty_uniform        $\leq 0.0181$     $\leq 0.0233$     1.0     $\leq 0.0221$     1.0     1

faulty_rare           $\leq 0.0019$     $\leq 0.0022$     1.0     $\leq 0.0022$     1.0     1

faulty_easy1          $\leq 0.0801$     $\leq 0.1007$     1.0     $\leq 0.0865$     1.0     1

faulty_ndecr          $\leq 0.0561$     $\leq 0.0723$     1.0     $\leq 0.0630$     1.0     1

faulty_walk           $\leq 0.0121$     $\leq 0.0173$     1.0     $\leq 0.0166$     1.0     1

legend cell align=left,
legend style=
fill opacity=0.8,
draw opacity=1,
text opacity=1,
anchor=south east,
log basis x=10,
log basis y=10,
tick align=outside,
tick pos=left,
x grid style=darkgray176,
xlabel style=align=center,
xlabel=Probability bound $p$
Program-Agnostic Neural Supermartingale,
xmin=0.001, xmax=1.2,
xtick style=color=black,
y grid style=darkgray176,
ylabel style=align=center,
ylabel=Probability bound $p$
Farkas' Lemma,
ymin=0.001, ymax=1.2,
ytick style=color=black,
[semithick, crimson2143940, mark=*, mark size=3, mark options=solid, only marks]
0.0842 1
0.0739 1
0.0553 1
0.0313 1
Farkas Fails
[semithick, crimson2143940, mark=x, mark size=3, mark options=solid, only marks]
0.0288 0.0991
0.0344 0.0991
0.0568 0.0991
0.0864 0.1819
0.1399 0.1819
0.0095 0.0091
0.0195 0.0181
0.0233 0.0181
0.0022 0.0019
0.1007 0.0801
0.0723 0.0561
0.0173 0.0121
0.0553 0.4546
Farkas Succeeds
[semithick, black, dashed]
0.0001 0.0001
1 1

Probability bounds generated using program-agnostic neural supermartingales and using Farkas' Lemma. The reference line $y=x$ allows one to see which approach outperforms the other and by how much: above the line means that neural supermartingales outperform Farkas' Lemma, below the line the opposite. Neural supermartingales can significantly outperform linear templates when a better bound exists, but otherwise achieve similar results. Our approach with program-aware neural supermartingales provides even better outcomes, compared with Farkas' Lemma.

The previous section develops a method for synthesising neural supermartingales. This section presents an empirical evaluation of the method, by testing it against a series of benchmarks. Each benchmark is tested ten times. A test is successful if a valid supermartingale is synthesised, and the proportion of successful tests is recorded. Further, the average bound from the valid supermartingales is also recorded, along with the average time taken by the learning and verification steps, respectively. This procedure is applied separately for program-aware and program-agnostic synthesis.
We use values $\beta_1 = 10, \beta_2 = 10^5, \beta_3 = 1$ in Eq. (<ref>) to define the priority ordering of the three terms in the loss function, which we find beneficial across our set of benchmarks.
To compare our method against existing work, we perform template-based synthesis of linear supermartingales using Farkas' Lemma. This requires deterministic invariants to overapproximate the reachable set of states, which may either be generated by abstract interpretation, or provided manually [48]. In our experiments, we provide a suitable invariant manually based on the guard of the loop, in some cases strengthening them with additional constraints by an educated guess.

It should be noted that our method is inherently stochastic. One reason is the random initialisation of the neural template's parameters in the learning phase. In program-agnostic synthesis, an additional source of randomness is the sampling of successor states. So that the results accurately reflect the performance of our method, the random seed for these sources of randomness is selected differently for each test. An additional source of non-determinism arises from the SMT solver Z3 as it generates counterexamples: this cannot be controlled externally. Benchmarks are run on a machine with an Nvidia A40 GPU, and
involve the assertion-violation analysis of programs created using the following two patterns:

Unreliable Hardware These are programs that execute on  unreliable hardware. The goal is to upper bound the probability that the program fails to terminate due to a hardware fault. A simple example is , whose source code is presented in <ref>, <Ref>, and which consists of a loop which may violate an assertion with small probability, modelling a hardware fault.

Robot Motion These programs model an agent (e.g., a robot) that moves within a physical environment. In these benchmarks, the uncertainty in control and sensing is modelled probabilistically. The environment contains a target region and a hazardous region. The goal is to upper bound the probability that the robot enters the hazardous region. We provide the program  as an example (<ref>, <ref>), which is a variant of repulse (<ref>) containing a modified assertion and initial state. The program models the motion of a robot in a one-dimensional environment, starting at ${\tt x} = 10$. The target region is where ${\tt x} < 0$, and hazardous region is where ${\tt x} > 100$. As with repulse, in each iteration there is an equal probability of ${\tt x}$ being decremented by 2, and $x$ being incremented by 1.

Several of the programs are based on benchmarks used in prior work focusing on other types of supermartingales [7, 15]. Additionally, there are several benchmarks that are entirely new.

Results showing the time taken in seconds to synthesise supermartingales by our method and Farkas' Lemma. For our method, we show the time taken during learning and verification.

0.88Benchmark     Farkas'     4cQuantitative Neural Certificates
Lemma     2cProgram-Agnostic     2cProgram-Aware
(lr)3-4 (lr)5-6
Learn Time     Verify Time     Learn Time     Verify Time

persist_2d            -          169.14     85.31     44.96     74.90

faulty_marbles        -                 114.24      29.23     15.86      28.68

faulty_unreliable     -                 123.85     45.48     18.34      33.97

faulty_regions        -                 17.92     35.85     17.55     32.38
cliff_crossing        0.11          134.61     19.02     21.27     29.07

repulse100                0.19      16.65      5.00      6.49      3.74

repulse100_uniform       0.19      21.28      14.18     -         -

repulse100_2d            0.12     122.92     64.54     15.75     47.70

faulty_varying        0.36      21.74      5.06      4.71      3.28

faulty_concave        0.39      49.12      13.37     13.49     7.82
fixed_loop            0.15      14.16      3.14      3.34      2.43

faulty_loop           0.16      25.52      3.81      3.73      2.66

faulty_uniform        0.34      20.20      1.91      6.75      1.33

faulty_rare           0.27      25.52      4.27      3.71      2.96

faulty_easy1          0.31      104.20     12.78     4.95      7.51

faulty_ndecr          0.33      104.89     9.06      5.37      4.66

faulty_walk           0.32      15.08      4.00      6.97      3.33

The results are reported in Table <ref>. The table is divided into three sections. The first section shows the benchmarks where Farkas' Lemma cannot be applied, and where only our method is capable of producing a bound. The second section shows examples where both methods are able to produce a bound, but our method produces a notably better bound. The third section shows benchmarks where both methods produce comparable bounds.

Dashes in the table indicate experiments where a valid supermartingale could not be obtained. In the case of Farkas' Lemma, there are several cases where there is no linear supermartingale for the benchmark. By contrast, program-agnostic and program-aware synthesis could be applied to almost all benchmarks. The one exception is  where only program-aware verification was unsuccessful: this is due to indicator functions in the post-expectation, which are not smooth and posed a problem for the optimiser. This benchmark underscores the value of program-agnostic synthesis, since it does not require embedding the explicit post-expectation in the loss function.

The first section of the benchmarks in <ref> demonstrates that our method produces useful results on programs that are out-of-scope for existing techniques. Furthermore, the success ratio of our method is high on all the benchmarks, which indicates its robustness. For the second and third sections (which consist of benchmarks to which Farkas' Lemma is applicable), the success ratio of our method is broadly maximal, which is to be expected, since these programs can also be solved by Farkas' Lemma.

In the second section of <ref>, we find more complex benchmarks where our method was able to significantly improve the bound from Farkas' Lemma. The smallest improvement was about 0.04, and the largest improvement was over 0.39. The intuition here is that neural templates allow more sophisticated supermartingales to be learnt, that can approximate how the reachability probability varies across the state space better than linear templates, and thereby yield tighter probability bounds.

The third section of <ref> consists of relatively simple benchmarks, where our method produces results that are marginally less tight in comparison to Farkas' Lemma. This is not surprising since our method uses neural networks consisting of a single neuron for these examples, owing to their simplicity. The expressive power of these networks is therefore similar to linear templates.

In summary, the results show that our method does significantly better on more complex examples, and marginally worse on very simple examples.
This is highlighted in <Ref>. Each point represents a benchmark. The position on the $x$-axis shows the probability bound obtained by our program-agnostic method, and the $y$-axis shows the probability bound obtained by Farkas' Lemma. Points above the line indicate benchmarks where neural supermartingales outperform Farkas' Lemma, and vice versa. The scale is logarithmic to emphasise order-of-magnitude differences.

Notice that in <ref> the program-aware algorithm usually yields better bounds than the program-agnostic algorithm, but the improvement is mostly marginal. This is in fact a strength of our method: our data-driven approach performs almost as well as one dependent on symbolic representations, which is promising in light of questions of scalability to more complex programs. We also include a breakdown of computation time (<Ref>)
which allows distinguishing between learning and verification overheads. Notably, Farkas' Lemma is significantly faster than our method, given that it relies on solving a convex optimisation problem via linear programming, whereas the synthesis of neural supermartingales is a non-convex optimisation problem that is addressed using gradient descent.

Having presented the experimental results, we shall further comment on some specific benchmarks.
The repulse100 program in <ref> is a variation of repulse presented in <ref>, but with the assertion changed to assert(x <= 100). While this is a small program, our method is still able to produce a significantly better result than Farkas' lemma, using a neural supermartingale with a single hidden layer consisting of three ReLU components that are summed together, which allows a convex piecewise linear function to be learnt.
The cliff_crossing program (<ref>, <ref>) is a further benchmark for which our method is capable of producing a significantly better bound. This is a 2 dimensional benchmark, for which we use a neural supermartingale that consists of two input neurons and four ReLU components, leading to a clear improvement compared to the linear supermartingale in the tightness of the probability bounds generated, as illustrated by <ref>.
Both repulse100 and cliff_crossing are benchmarks that use neural supermartingales with a single hidden layer. An example that uses two hidden layers is faulty_concave (<ref>, <ref>), in which there are two distinct regions of the state space, one of which has a significantly higher reachability probability than the other. We find that neither a linear template nor a single-layer neural supermartingale is able to exploit this conditional behaviour, each of which yield an overly conservative certificate, but that a neural supermartingale with two hidden layers is able to more tightly approximate the reachability probability in each of the two regions, which is illustrated in <ref>.
Further discussion and the source code of these case studies are presented in <ref>.

[mesh/ordering=y varies,colormap/greenyellow,width=2.9in]
3 [surf,mesh/rows=11] file graphs/cliff_crossing_linear.dat;

*Linear supermartingale.

[mesh/ordering=y varies,colormap/cool,width=2.9in]
3 [surf,mesh/rows=11] file graphs/cliff_crossing_neural.dat;

*Neural supermartingale.

Supermartingale functions for the  benchmark as generated using Farkas' Lemma (on the left) and using neural supermartingales (on the right). The right hand figure illustrates the tighter bounds obtainable through the use of neural templates.

§ RELATED WORK

The formal verification of probabilistic programs using supermartingales is a well-studied topic. Early approaches to introduce this technique applied them to almost-sure termination analysis of probabilistic programs [15], which allowed several extensions to polynomial programs, programs with non-determinism,
lexicographic and modular termination arguments, and
persistence properties [6, 17, 18, 21, 25, 31]. All these methods relied on symbolic reasoning algorithms for synthesising supermartingales, that leveraged theories based on Farkas' lemma for the synthesis of linear certificates, and Putinar's Positivestellensatz and sum-of-square methods for the synthesis of polynomial certificates. While these methods are the state-of-the-art for many existing problem instances in literature, to achieve the strong guarantees that they provide (such as completeness for the specific class of programs they target), they must necessarily introduce restrictions on the class of programs to which they are applicable, and the form of certificates that they derive.
Moreover, symbolic methods need externally provided invariants that are stronger than $\mathbb{R}^n$ to enforce non-negativity in the case of linear certificates, as we illustrate in <ref>. Also, symbolic methods for the synthesis of polynomial certificates require compact deterministic invariants to operate.

The use of neural networks to represent certificates has allowed many of these restrictions to be lifted. In the context of the analysis of probabilistic programs, neural networks were first applied to certify positive almost-sure termination [3]. This approach lent itself to a wider range of formal verification questions for stochastic dynamical models, from stability and safety analysis to controller synthesis [36, 20, 37]. These data-driven inductive synthesis techniques for supermartingales have also been extended to machine learning techniques other then deep learning, such as piecewise linear regression and decision tree learning [13, 12].

In this paper, we further extend data-driven synthesis of neural supermartingale certificates to quantitative verification questions. The correctness of our approach builds upon the theory of non-negative repulsing supermartingales [48], as formulated in <ref>.
Our experiments have demonstrated that neural certificates attain comparable results on programs that are amenable to symbolic analysis (such as those in the third section of <ref>), while surpassing symbolic methods on more complex programs that are either out-of-scope or yield overly conservative bounds when existing techniques are applied (such as those in the first and second section of <ref>).

§ CONCLUSION

We have presented a data-driven framework for the quantitative verification of probabilistic models
that leverage neural networks to represent supermartingale certificates.
Our experiments have shown that neural certificates are applicable to a wider range of probabilistic models than was previously possible using purely symbolic techniques. We also illustrate that on existing models our method yields certificates of better or comparable quality than those produced by symbolic techniques for the synthesis of linear supermartingales. This builds upon the ability of neural networks to approximate non-linear functions, while satisfying the constraints imposed by <ref> without the need for supporting deterministic invariants to be provided externally. Our method applies to quantitative termination and assertion-violation analysis for probabilistic programs, as well as safety and invariant verification for stochastic dynamic models. We imagine extensions to further quantitative verification questions, such as temporal properties beyond reachability [11], and bounding expected accrued costs [42, 52, 51].

§ PROOFS

§.§ Notation

We denote the set of non-negative reals as $\mathbb{R}_{\geq 0}$ and the set of non-negative extended reals as $\mathbb{R}_{\geq 0} \cup \{ \infty \}$.
Given a function $f : X \to Y$, we define the image $f[A]$ of a set $A \subseteq X$ under $f$ as the set $f[A]= \{ f(x) ~:~ x \in A \}$;
we define the pre-image $f^{-1}[B]$ of a set $B \subseteq Y$ under $f$ by $f^{-1}[B] = \{x \in X : f(x) \in B \}$.

§.§ Preliminaries of Measure Theory

Let $(X, \mathcal{A})$ be a measurable space.
We say that a sequence $\{ f_n \}_{n \in \mathbb{N}}$
of non-negative $(\mathcal{A}, \mathcal{B}(\mathbb{R}_{\geq 0} \cup \{ \infty\}))$-measurable functions $f_n : X \to \mathbb{R}_{\geq 0}  \cup \{ \infty \}$  converges monotonically to $f$ if the following two conditions hold:

* the sequence converges pointwise to $f$, i.e,
$\lim_{n \to \infty} f_n(x) = f(x)$ for all $x \in X$, and

* the convergence is monotone, i.e.,
$f_{n}(x) \leq f_{n+1}(x)$ for all $n \in \mathbb{N}$ and $x \in X$.

Let $(X, \mathcal{A}, \mu)$ be a measure space, and $A \in \mathcal{A}$ be a measurable set, then we note that

\begin{equation}
\mu(A) = \int_X \mathbf{1}_{A}(x)~\mu(\mathrm{d}x) \label{eq:int}
\end{equation}

namely, a set's measure is equal to the integral of its indicator function with respect to that measure.
Furthermore, integrals satisfy the following theorem.

Let $(X, \mathcal{A}, \mu)$ be a measure space.
Let $\{ f_n \}_{n \in \mathbb{N}} $ be a sequence of non-negative $(\mathcal{A}, \mathcal{B}(\mathbb{R}_{\geq 0} \cup \{ \infty\}))$-measurable functions $f_n : X \to \mathbb{R}_{\geq 0}  \cup \{ \infty \}$ that converge monotonically to $f$.

\begin{equation*}
\int_{X} ~f(x) ~\mu(\mathrm{d}x) = \lim_{n \to \infty} \int_{X} ~f_n(x) ~\mu(\mathrm{d}x).
\end{equation*}

Let $(X, \mathcal{A}, \mu)$ be a measure space.
Let $\{ A_n \}_{n \in \mathbb{N}}$ be a sequence of measurable subsets $A_i \in \mathcal{A}$ such that $A_n \subseteq A_{n+1}$ for all $n \in \mathbb{N}$. Let $A = \cup \{ A_n : n \in \mathbb{N} \}$. Then

\begin{equation*}
\mu( A) =
\lim_{n \to \infty} \mu(A_n)
\end{equation*}

It follows from Theorem <ref> applied to the indicator functions of the sets $A_n$ and the indicator function of $A$. We begin by showing that the sequence of indicator functions
$\{ \mathbf{1}_{A_n} \}_{n \in \mathbb{N}}$ converges monotonically to $\mathbf{1}_{A}$.

* First, to show that the sequence converges pointwise to $\mathbf{1}_{A}$, we consider the two cases where $x \in A$, and where $x \not\in A$. In the case that $x \in A$, then there exists some $N \in \mathbb{N}$ for which for $x \in A_N$, since $A = \cup \{ A_n : n \in \mathbb{N} \}$. Then, since $A_n \subseteq A_{n + 1}$, we also have that $x \in A_n$ for all $n \geq N$. This implies that $\lim_{n \to \infty} \mathbf{1}_{A_n}(x) = \mathbf{1}_{A}(x)$. In the case that $x \notin A$, then none of the sets $A_n$ contain $x$.
This implies that $\lim_{n\to\infty}\mathbf{1}_{A_n}(x)$ and $\mathbf{1}_{A}(x)$ are both equal to zero.

* Second, to show that the convergence is monotone, we simply observe that   $A_n \subseteq A_{n + 1}$
implies that $\mathbf{1}_{A_{n}}(x) \leq \mathbf{1}_{A_{n+1}}(x)$.

Finally, we invoke the Monotone Convergence Theorem as follows:

\begin{align*}
\mu(A) &= \int_X \mathbf{1}_{A}(x) ~\mu(\mathrm{d}x) &\text{By ~\eqref{eq:int}}\\
&= \lim_{n \to \infty} \int_X \mathbf{1}_{A_n}(x) ~\mu(\mathrm{d}x) &\text{By \cref{thm:mct}} \\
&= \lim_{n \to \infty} \mu( A_n ) &\text{By~\eqref{eq:int}} &\qedhere
\end{align*}

§.§ Preliminaries of Domain Theory

Let $(P, \sqsubseteq)$ be a partially ordered set. We say that $(P, \sqsubseteq)$ is a complete partially ordered set (CPO) if
every increasing chain $p_1 \sqsubseteq p_2 \sqsubseteq \dots $ has a least upper bound $\sqcup \{ p_i ~:~ i \in \mathbb{N} \}$ in $P$.
We say that a CPO is pointed if it has a bottom element $\bot$.

We say that an endofunction $F \colon P \to P$ on a CPO $(P, \sqsubseteq)$ is
Scott-continuous if it is

* monotonic, i.e., for all $p \sqsubseteq p' \in P$, we have that $F(p) \sqsubseteq F\left(p^{\prime}\right)$,  and

* preserves least upper bounds of countable increasing chains, i.e., for all chains $p_0 \sqsubseteq p_1 \sqsubseteq \dots $ in $P$, we have
\sup~\{
F(p_i) ~:~ i \in \mathbb{N}
\}
p_i ~:~ i \in \mathbb{N}
\}).$

For a function $F : P \to P$ we define, for $x \in P$, $F^0(x) = x$ and $F^{n+1}(x) = F(F^{n}(x))$.
A pre-fixpoint of $F$ is an element $p$ of $P$ such that $F(p) \sqsubseteq p$; a fixpoint of $F$ is an element $f$ of $P$ such that $F(f) = f$; a least fixpoint is a fixpoint that is least w.r.t. $\sqsubseteq$.

Let $F : P \to P$ be a Scott-continuous endofunction on a pointed CPO. Then, $F$ has a least fixpoint, which is $\sup~\{ F^n(\bot) ~:~ n \in \mathbb{N} \}$.

§.§ Proof of Lemma <ref>

For every $i$ we have that the pre-image of
any measurable set $R  \in \mathcal{B}(\mathbb{R}^n)$ under $X_i^{s} : (\Omega, \mathcal{F}) \to (\mathbb{R}^n, \mathcal{B}(\mathbb{R}^n))$, denoted $(X_{i}^s)^{-1}[R]$, is measurable since $X_i^{s}$ is a measurable function <cit.>. It follows that $\Reach_t^{s}(A)$ is measurable because it can be expressed as the following finite intersection:

\begin{equation*}
\Reach_t^{s}(A) =
\cap\{ (X_i^{s})^{-1}[\mathbb{R}^n \setminus A] ~\colon~ i = 0, \dots, t-1\}
\cap (X_t^s)^{-1} [A].
\end{equation*}

It is a standard result (see <cit.> and [5]) that the probability measure of the event $\Reach_{{t}}^{s_0}(A)$ with $t > 0$ can be expressed as an integral with respect to the transition kernel of the
Markov chain:

\begin{align*}
\mathbb{P}[\Reach_{{t}}^{s_0}(A)] &= \mathbf{1}_{\mathbb{R}^n \setminus A}(s_0)
\int_{s_1 \in \mathbb{R}^n \setminus A} \cdots
\int_{s_{t-1} \in \mathbb{R}^n \setminus A}  \int_{s_{t} \in A}
~ \prod_{i=1}^t T(s_{i-1}, \mathrm{d}s_{i})\\
&= \int_{s_1 \in \mathbb{R}^n} \cdots \int_{s_{t} \in \mathbb{R}^n}
\left( \prod_{i=0}^{t-1} \mathbf{1}_{\mathbb{R}^n \setminus A}(s_i) \right) \mathbf{1}_{A}(s_t) \prod_{i=1}^{t} T(s_{i-1}, \mathrm{d}s_{i})
\end{align*}

As a consequence, we can express the probability measure of
$\Reach_{t+1}^{s_0}(A)$ as follows:

\begin{multline*}
\mathbb{P}[\Reach_{t+1}^{s_0}(A)] = \mathbf{1}_{\mathbb{R}^n \setminus A}(s_{0})
\int_{s_1 \in \Real^n}
\\
\underbrace{\int_{s_2 \in \Real^n} \cdots \int_{s_{t+1} \in \Real^n}
\left( \prod_{i=1}^t \mathbf{1}_{\mathbb{R}^n \setminus A}(s_i) \right)
\mathbf{1}_{A}(s_{t+1})
\left(\prod_{i=2}^{ t+1 } T(s_{i-1}, \mathrm{d}s_{i}) \right)}_{\mathbb{P}[\Reach_{t}^{s_1}(A)]}
T(s_0, \mathrm{d}s_1)
\end{multline*}

Then, after substituting $s_0$ with $s$ and $s_1$ with $s'$, we derive Eq. (<ref>) as follows:

\begin{align*}
\mathbb{P}[\Reach_{{t+1}}^{s}(A)]
\mathbf{1}_{\mathbb{R}^n \setminus A}(s)
\cdot
\int_{\Real^n}
\mathbb{P}[\Reach_{{t}}^{s'}(A)] ~T(s, \mathrm{d}s')
& \\
\mathbf{1}_{\mathbb{R}^n \setminus A}(s)
\cdot
\mathbb{X}[\lambda s' \mathpunct{.}
\mathbb{P}[\Reach_{{t}}^{s'}(A)]
&\text{By \eqref{eqn:postex}}
\end{align*}

Finally, Eq. (<ref>), which states that $\mathbb{P}[\Reach_{0}^{s}(A)] = \mathbf{1}_A(s)$ directly follows
from the definition of Reach, as
$\Reach_{0}^{s}(A) = \Omega$ when $s \in A$ and
$\Reach_{0}^{s}(A) = \emptyset$ when $s \not\in A$.

§.§ Proof of Lemma <ref>

Since the sets in the family $\{ \Reach^{s}_{t} \colon t \in \Nat\}$ are disjoint, we note that the union in (<ref>) that defines $\Reach^{s}_{\leq t}$ is a union of disjoint sets. It follows that, by countable additivity of measure and disjointness of the union in (<ref>), its probability measure is given by

\begin{equation}\label{eqn:reach-sum}
\mathbb{P}[\Reach_{\leq t}^{s}(A)] =
\sum_{i=0}^t  \mathbb{P}[\Reach_{t}^{s}(A)]
\end{equation}

Consequently and thanks to <ref>, we have

\begin{align*}
\mathbb{P}[\Reach_{\leq {t+1}}^{s}(A)]
&= \mathbb{P}[\Reach_{{0}}^{s}(A)] + \sum_{i = 1}^{t + 1} \mathbb{P}[\Reach_{{i}}^{s}(A)]
&= \mathbf{1}_{A}(s) + \sum_{i = 1}^{t + 1} \mathbb{P}[\Reach_{{i}}^{s}(A)] &\text{By~\eqref{eqn:standard-result-zero}}\\
&= \mathbf{1}_{A}(s) + \sum_{i = 0}^{t} \mathbb{P}[\Reach_{{i + 1}}^{s}(A)]\\
\mathbf{1}_{A}(s) + \sum_{i = 0}^{t}
\mathbf{1}_{\mathbb{R}^n \setminus A}(s) \cdot
\mathbb{X}[\lambda s' \mathpunct{.}
\mathbb{P}[\Reach_{{i}}^{s'}(A)]](s)
&\text{By \eqref{eqn:standard-result-reach}}\\
\mathbf{1}_{A}(s) + \sum_{i = 0}^{t}\left(
\mathbf{1}_{\mathbb{R}^n \setminus A}(s) \cdot
\int_{\Real^n}
\mathbb{P}[\Reach_{{i}}^{ s^\prime }(A)] ~T(s, \mathrm{d}s')
\right) &\text{By \eqref{eqn:postex}}\\
\mathbf{1}_{A}(s) +
\mathbf{1}_{\mathbb{R}^n \setminus A}(s) \cdot
\int_{\Real^n}
\left( \sum_{i = 0}^{t} \mathbb{P}[\Reach_{{i}}^{s'}(A)] \right) ~T(s, \mathrm{d}s')
\\
\mathbf{1}_{A}(s) +
\mathbf{1}_{\mathbb{R}^n \setminus A}(s) \cdot
\int_{\Real^n}
\mathbb{P}[\Reach_{{\leq t}}^{s'}(A)] ~T(s, \mathrm{d}s')
\mathbf{1}_{\Real^n \setminus A}(s)
\cdot
\mathbb{X}[\lambda s' \mathpunct{.}
\mathbb{P}[\Reach_{{\leq t}}^{s'}(A)]](s)
\end{align*}
Furthermore, we also have that

\begin{align*}
\mathbb{P}[\Reach_{\leq 0}^{s_0}(A)]
&= \mathbb{P}[\Reach_{0}^{s_0}(A)]
= \mathbf{1}_{A}(s_0)
&\text{By~\eqref{eqn:reach-sum} and \eqref{eqn:standard-result-zero}}.
\end{align*}

§.§ Proof of Lemma <ref>

We recall from Eq. (<ref>) that $\Reach_{\leq t}^{s}(A) = \cup \left\{\Reach_{i}^{s}(A) \colon i \leq t \right\}$,
and we note that this implies that
$\Reach_{\leq t}^{s}(A) \subseteq \Reach_{\leq t+1}^{s}(A)$.
Then, we recall from Eq. (<ref>) that $\Reach_{\mathrm{fin}}^s(A) = \cup \left\{\Reach_{\leq t}^{s}(A) \colon t \in \Nat\right\}$. Finally, we apply <ref> to conclude $\mathbb P [\Reach_{\mathrm{fin}}^s(A)] = \lim_{t \to \infty}\mathbb P [\Reach_{\leq t}^s(A)]$.

§.§ Proof of Theorem <ref>

We adapt the soundness proof of non-negative repulsing supermartingales [48] to neural certificates (which do not require an externally-provided deterministic invariant to restrict their domain) in our setting of probabilistic systems (i.e. we do not consider non-deterministic transition dynamics).

Let $\mathcal{V} \subset \Real^n \to \Real_{\geq 0} \cup \{ \infty \}$ be the set of $(\mathcal{B}(\mathbb{R}^n), \mathcal{B}(\mathbb{R}_{\geq 0} \cup \{ \infty\}))$-measurable functions.
Let the ordering relation $\leq$ between functions be defined pointwise as

\begin{equation}\label{eqn:order}
U \leq V \Leftrightarrow \forall s \in \Real^n \mathpunct{.} U(s) \leq V(s)
\end{equation}

Then $(\mathcal{V}, \leq)$ is a pointed CPO and $0$ (defined as the constant zero function) is its bottom element. Furthermore,
for every countable increasing chain $\{V_i \}_{i \in \mathbb{N}}$, the following holds:

\begin{equation}\label{eqn:lublim}
(\sup~ \{V_i \colon i \in \Nat \})(s)
= \sup~\{ V_i(s) : i \in \Nat \}
= \lim_{i \to \infty} V_i(s).
\end{equation}

The fact that the relation $\leq$ on functions is reflexive, transitive and anti-symmetric follows from the properties of the natural ordering relation $\leq$ on the set $\mathbb{R}_{\geq 0} \cup \{ \infty \}$. This shows that $(\mathcal{V}, \leq)$ is a partial order. The constant function $0$ is the bottom element of the partial order because $\mathcal{V}$ contains only non-negative functions.

To show that the partial order is also complete, we need to show that the least upper bound of every countable increasing chain $\{V_i \colon i \in \mathbb{N}\}$ in $(\mathcal{V}, \leq)$ is in $\mathcal{V}$;
in other words, the least upper bound is $(\mathcal{B}(\mathbb{R}^n), \mathcal{B}(\mathbb{R}_{\geq 0} \cup \{ \infty\}))$-measurable.
To this goal, we first show that
$(\sup~ \{V_i \colon i \in \Nat \})(s) = \sup~\{ V_i(s) : i \in \Nat \}$, and then that this is measurable.
We observe that any function $V_{\sup}$ that maps $s \in \Real^n$ to $\sup~\{ V_i(s) : i \in \Nat \}$ must satisfy the following two properties:

* $\forall s \in \Real^n . \forall i \in \mathbb{N} . V_i(s) \leq V_{\sup}(s)$

* $\forall s \in \Real^n . \forall b \in \mathbb{R}_{\geq 0} \cup \{ \infty\} .~ [\forall i \in \mathbb{N} .~  V_i(s) \leq b] \implies V_{\sup}(s) \leq b$.

Notably, the second statement applies to all functions whose co-domain has type $\Real_{\geq 0} \cup \{ \infty \}$, including those with domain $\Real^n$. Therefore

* $\forall s \in \Real^n .
\forall U \in \Real^n \to \Real_{\geq 0} \cup \{ \infty \} .~ [
\forall i \in \mathbb{N} .~  V_i(s) \leq U(s)] \implies
V_{\sup}(s) \leq U(s)$.

Since the set of measurable functions $\mathcal{V}$ is a subset of $\Real^n \to \Real_{\geq 0} \cup \{ \infty \}$, this implies

* $\forall s \in \Real^n .\forall U \in \mathcal{V} .~ [
\forall i \in \mathbb{N} .~  V_i(s) \leq U(s)] \implies
V_{\sup}(s) \leq U(s)$.

Finally, we note that this implies the following sentence:

* $\forall U \in \mathcal{V} .~ [
\forall s \in \Real^n .
\forall i \in \mathbb{N} .
V_i(s) \leq U(s)] \implies
\forall s \in \Real^n .
V_{\sup}(s) \leq U(s)$.

Therefore, we have that function $V_{\sup}$ also satisfies the following two properties:

* $\forall i \in \mathbb{N} . \forall s \in \Real^n  . V_i(s) \leq V_{\sup}(s)$

* $\forall U \in \mathcal{V} .~ [
\forall i \in \mathbb{N} .
\forall s \in \Real^n .
V_i(s) \leq U(s)] \implies
\forall s \in \Real^n .
V_{\sup}(s) \leq U(s)$.

Using (<ref>), we can rewrite these two sentences in terms of the ordering relation over functions:

* $\forall i \in \mathbb{N} . V_i \leq V_{\sup}$

* $\forall U \in \mathcal{V} .~ [
\forall i \in \mathbb{N} .~  V_i \leq U] \implies
V_{\sup} \leq U$.

But these are precisely the statements that function $\sup~\{V_i \colon i \in \mathbb{N} \}$ must satisfy. Since the supremum must be unique, this implies that
$V_{\sup}$ and $\sup~\{V_i \colon i \in \mathbb{N} \}$ must be equal.
Then, to show that $V_{\sup}$ is measurable, we rely on the standard result that the pointwise supremum of any countable family of measurable functions is also measurable <cit.>.

Finally, since $\{V_i(s)\}_{i \in \mathbb{N}}$ is a non-decreasing sequence in the extended non-negative reals, we have by the monotone convergence property of non-decreasing real sequences <cit.> that

\begin{equation*}
\lim_{i \to \infty} V_i(s) = \sup~\{ V_i(s) \colon i \in \mathbb{N}\}.\qedhere
\end{equation*}

The post-expectation operator $\mathbb{X}$ is Scott-continuous w.r.t. $(\mathcal{V}, \leq)$.

We show that $\mathbb{X}$ (defined in Eq. (<ref>)) satisfies the following two properties.

* $\mathbb{X}$ is monotonic. Suppose $U,V \in \mathcal{V}$ with $U \leq V$. Then for all $s \in \mathbb{R}^n$:

\begin{align*}
\mathbb{X}[U](s) &= \int_{\mathbb{R}^n} U(s^\prime) ~T(s, \mathrm{d}s^\prime)\\
\int_{\mathbb{R}^n}
V(s^\prime) ~T(s, \mathrm{d}s^\prime)
&&\text{By $U \leq V$}\\
&= \mathbb{X}[V](s).
\end{align*}

* $\mathbb{X}$ preserves least upper bounds of countable increasing chains. Let $V_0 \leq V_1 \leq \ldots $ be a countable increasing chain in $\mathcal{V}$. As $\mathbb{X}$ is monotone, $\mathbb{X}[V_0] \leq \mathbb{X}[V_1] \leq \ldots$ is also a countable increasing chain.
Then for all $s \in \mathbb{R}^n$:

\begin{align*}
(\sup~ \{
\mathbb{X}[V_i] \colon i \in \mathbb{N}
\})(s)
&= \lim_{i \to \infty} \mathbb{X}[V_i](s)
&\text{By~\eqref{eqn:lublim}, \cref{lem:cpo}}\\
&= \lim_{i \to \infty}
\int_{\Real^n} V_i(s^\prime)~T(s, \mathrm{d}s^\prime)
&\text{By \eqref{eqn:postex}}
\\
&= \int_{\Real^n} (\sup~\{ V_i \colon i \in \mathbb{N} \})(s^\prime)~T(s, \mathrm{d}s^\prime)
&\text{By \cref{thm:mct}}
\\
&= \mathbb{X}[\sup~\{ V_i \colon i \in \mathbb{N} \}](s)
\end{align*}

This involves an application of the Monotone Convergence Theorem (<ref>), as by <ref> we have that $\{V_i\}_{i \in \mathbb{N}}$ is a sequence of non-negative measurable functions that converge monotonically to the function $\sup~\{ V_i \colon i \in \mathbb{N}\}$.

Let $\Phi : \mathcal{V} \to \mathcal{V}$ be defined as the higher-order function

\begin{equation}\label{def:phi}
\Phi(V)(s) =
\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[V](s)
\end{equation}

Then, $\Phi$ is Scott-continuous w.r.t. $(\mathcal{V}, \leq)$.

We show that

* $\Phi$ is monotonic.
Let $s \in \Real^n$ be arbitrary and assume $U \leq V$:

\begin{align*}
\Phi(U)(s) &=
\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[U](s)\\
\leq \mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[V](s) &\text{By monotonicity of $\mathbb{X}$, \cref{lem:nexttime-cont}}\\
&= \Phi(V)(s),
\end{align*}

therefore $\Phi(U) \leq \Phi(V)$.

* $\Phi$ preserves least upper bounds of countable increasing chains. To show that $\Phi$ preserves least upper bounds of countable increasing chains, suppose $\{ V_i \colon i \in \mathbb{N}\}$ is a countable increasing chain in $(\mathcal{V}, \leq)$.
Let $s \in \Real^n$ be arbitrary:

\begin{align*}
\Phi(\sup~ \{V_i \colon i \in \Nat \})(s)
\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[
\sup~ \{V_i \colon i \in \Nat \}
](s) &\text{By \eqref{def:phi}}\\
\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot
(\sup~ \{\mathbb{X}[V_i] \colon i \in \Nat \}) (s) &\text{By \cref{lem:nexttime-cont}}\\
\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot
(\sup~\{\mathbb{X}[V_i](s) \colon i \in \Nat \}) &\text{By \eqref{eqn:lublim}}\\
\sup~ \{
(\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s)
\cdot \mathbb{X}[V_i](s)) \colon i \in \Nat
\}\\
\sup~ \{
\Phi[V_i](s) \colon i \in \Nat
\} &\text{By \eqref{def:phi}}\\
&= \left(\sup~ \{ \Phi(V_i) \colon i \in \Nat \}\right)(s)
&\text{By \eqref{eqn:lublim},}
\end{align*}

therefore $\Phi(\sup~ \{ V_i \colon i \in \mathbb{N}\}) = \sup~ \{ \Phi(V_i) \colon i \in \mathbb{N} \}$.

$\lambda s \mathpunct{.} \mathbb{P}[\Reach_{\mathrm{fin}}^{s}(A)]$ is the least fixpoint of $\Phi$.

Let $\Phi^t$ denote the $t$-th application of $\Phi$ to its argument.
By <ref>, $\Phi$ is a Scott-continuous endofunction on $(\mathcal{V}, \leq)$. Then, as a consequence of <ref>, $\Phi$ has a least fixpoint and this is equal to $\sup~ \{ \Phi^n(0) \colon n \in \mathbb{N}\}$.
Therefore, we need to prove

\begin{equation}
\lambda s . \mathbb{P}[\Reach_{\mathrm{fin}}^{s}(A)]
= (\sup~ \{ \Phi^n(0) \colon n \in \mathbb{N}\})
\label{eqn:reachlublambda}
\end{equation}

To prove Eq. (<ref>) we show first that for all $t \in \Nat$ and $s \in \Real^n$, we have

\begin{equation}
\mathbb{P}[\Reach_{\leq t}^{s}(A)]
= \Phi^{t+1}(0)(s)\label{eqn:kleene-1}
\end{equation}

We proceed by induction on $t \in \Nat$, that is, we let $s \in \mathbb{R}^n$ and show

* $\mathbb{P}[\Reach_{\leq 0}^{s}(A)]
= \Phi^{1}(0)(s)$, and

* $\mathbb{P}[\Reach_{\leq t}^{s}(A)]
= \Phi^{t+1}(0)(s)$
$\mathbb{P}[\Reach_{\leq t+1}^{s}(A)]
= \Phi^{t+2}(0)(s)$.

For the base case, we have

\begin{align*}
\mathbb{P}[\Reach_{\leq 0}^{s}(A)]
&= \mathbf{1}_A(s) &\text{By \eqref{eqn:reachbeforezero}}\\
&= \mathbf{1}_A(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[0](s)
&\text{Because $\mathbb{X}[0] = 0$}\\
&= \Phi^1(0)(s) &\text{By \eqref{def:phi}.}
\end{align*}

For the inductive case, we leverage <ref> as follows:

\begin{align*}
\mathbb{P}[\Reach_{\leq t + 1}^{s}(A)]
\mathbf{1}_A(s)
\mathbf{1}_{\Real^n \setminus A}(s)
\cdot
\mathbb{X}[
\lambda s^\prime .
\mathbb{P}[\Reach_{\leq t}^{s^\prime}(A)]
&\text{By \eqref{eqn:reachafterzero}}    \\
&= \Phi(\lambda s^\prime .  \mathbb{P}[\Reach_{\leq t}^{s^\prime}(A)])(s) &\text{By \eqref{def:phi}}\\
&= \Phi(\Phi^{t+1}(0))(s)&\text{By IH}\\
&= \Phi^{t+2}(0)(s)
\end{align*}

Next, we show that (<ref>) is a consequence of <ref>, noting that $\{\Phi^t(0) \}_{ t \in \mathbb{N}}$ is a countable increasing chain in  $(\mathcal{V},\leq)$ and $\Phi$ is monotonic:

\begin{align*}
\mathbb{P}[\Reach_{\mathrm{fin}}^{s}(A)] &= \lim_{t \to \infty} \mathbb{P}[\Reach_{\leq t}^{s}(A)] &\text{By \eqref{eqn:reachfinlim}}\\
&= \lim_{t \to \infty} \Phi^{t + 1}(0)(s) &\text{By \eqref{eqn:kleene-1}}\\
&= \lim_{t \to \infty} \Phi^{t}(0)(s) \\
&= (\sup~\{ \Phi^{t}(0) \colon t \in \Nat \})(s) &\text{By~\eqref{eqn:lublim}, \cref{lem:cpo}} &\qedhere
\end{align*}

Every function $V \in \mathcal{V}$ that satisfies  (<ref>) and (<ref>) is a pre-fixpoint of

Let us recall that $\Phi(V)(s) = \mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[V](s)$ by definition (<ref>). Then, we let $s \in \Real^n$ be arbitrary,
and show that $\Phi(V)(s) \leq V(s)$. We proceed by cases:

\begin{align*}
&\text{Case $s \in A$: }~\Phi(V)(s)=
\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[V](s) = 1 \leq V(s) &\text{By \eqref{eq:spec_ind}}\\
&\text{Case $s \notin A$: }~
\Phi(V)(s)=
\mathbf{1}_{A}(s) + \mathbf{1}_{\Real^n \setminus A}(s) \cdot \mathbb{X}[V](s)
= \mathbb{X}[V](s)
\leq V(s) &\text{By \eqref{eq:spec-dec}}&\qedhere
\end{align*}

Finally, to conclude <ref>, we suppose that $V \in \mathcal{V}$ satisfies (<ref>) and (<ref>).
We recall from <ref> that
$\lambda s \mathpunct{.} \mathbb{P}[\Reach_{\mathrm{fin}}^{s}(A)]$ is the least fixpoint of $\Phi$
and, equivalently, the least pre-fixpoint of $\Phi$.
Since from <ref> we have that $V$ is a pre-fixpoint of $\Phi$, then $\lambda s \mathpunct{.} \mathbb{P}[\Reach_{\mathrm{fin}}^{s}(A)] \leq V$ holds true,
which corresponds to the statement of <ref>.

§ CASE STUDIES

In this section, we elaborate on the source code of the case studies discussed in <ref>:

* The program  (<ref>) simulates a program executing on unreliable hardware. A variable  is initialised to 0. In each successful iteration, it is incremented by 1 with probability 0.5. A hardware fault is modelled with the  flag being set to 1.

* repulse100 (<ref>): models the motion of a robot in a one-dimensional environment, starting at position ${\tt x} =10$ (<ref>).

*  (<ref>) is a benchmark for which our method is capable of producing a significantly better bound that that of Farkas' Lemma. This program models a robot moving along a road that is adjacent to a cliff. The $x$-coordinate models how far along the road the robot is, with $x = 0$ being the beginning, and $x = 50$ being the end. Then, the $y$-coordinate models how close to the cliff the robot is, with $y = 0$ being the furthest point away from the cliff, and $y = 10$ being the closest point to the cliff still on the road. Our method computes an upper bound on the probability of the assertion assert(y <= 10) being violated prior to the robot reaching ${\tt x > 50}$.

*  (<ref>): the variable i is initialised to 1. In each successful iteration, it is incremented with probability 0.5 until 10 is reached. In each iteration, there is a small probability of a hardware fault. Unlike , the probability of failure changes depending on the value of $i$. When $1 \leq i \leq 3$, the probability of failure is 0.001. However, when $4 \leq i \leq 9$, the probability of failure is 0.01, an order of magnitude higher. This is an example for which a neural supermartingale with at least two hidden layers is required, and for which neither a linear supermartingale nor a neural supermartingale with a single hidden layer suffice.

[caption=faulty_loop, label=lst:faulty]
error = 0;
i = 0;
while i < 10 do
assert(error == 0);
p   Bernoulli(0.999);
if p == 1 then
q   Bernoulli(0.5);
if q == 1 then
i += 1
error = 1

[caption=repulse100 , label=lst:repulse]
x = 10;
while x >= 0 do
assert(x <= 100);
p   Bernoulli(0.5);
if p == 1 then
x -= 2
x += 1

[caption=cliff_crossing ,label=lst:cliff]
x = 0;
y = 0;
while x <= 50 do
assert(y <= 10);
p   Bernoulli(0.9);
x += 1;
if p == 0 then
y += 1
if y >= 1 then
y -= 1

[caption=faulty_concave ,label=lst:concave]
i = 1;
error = 0;
while i < 10 do
assert(error == 0);
if i < 4 then
p   Bernoulli(0.999)
p   Bernoulli(0.99)
if p == 1 then
q   Bernoulli(0.5);
if q == 1 then
i += 1
error = 1

Linear and neural supermartingale certificates generated using Farkas' Lemma, and program-agnostic neural synthesis, respectively, for the   and  benchmarks. The values of the certificates at the initial state (dashed line) shows that tighter upper bounds are generated using neural certificates compared to Farkas' Lemma.

[1]
Alessandro Abate, Daniele Ahmed, Mirco Giacobbe, and Andrea Peruffo.
Formal synthesis of Lyapunov neural networks.
IEEE Control. Syst. Lett., 5(3):773–778, 2021.

[2]
Alessandro Abate, Alec Edwards, Mirco Giacobbe, Hashan Punchihewa, and Diptarko Roy.
Quantitative verification with neural networks.
In CONCUR, volume 279 of LIPIcs, pages 22:1–22:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2023.

[3]
Alessandro Abate, Mirco Giacobbe, and Diptarko Roy.
Learning probabilistic termination proofs.
In CAV (2), volume 12760 of Lecture Notes in Computer Science, pages 3–26. Springer, 2021.

[4]
Alessandro Abate, Joost-Pieter Katoen, and Alexandru Mereacre.
Quantitative automata model checking of autonomous stochastic hybrid systems.
In HSCC, pages 83–92. ACM, 2011.

[5]
Alessandro Abate, Maria Prandini, John Lygeros, and Shankar Sastry.
Probabilistic reachability and safety for controlled discrete time stochastic hybrid systems.
Automatica, 44(11):2724 – 2734, 2008.

[6]
Sheshansh Agrawal, Krishnendu Chatterjee, and Petr Novotný.
Lexicographic ranking supermartingales: an efficient approach to termination of probabilistic programs.
Proc. ACM Program. Lang., 2(POPL):34:1–34:32, 2018.

[7]
Christophe Alias, Alain Darte, Paul Feautrier, and Laure Gonnord.
Multi-dimensional rankings, program termination, and complexity bounds of flowchart programs.
In SAS, volume 6337 of Lecture Notes in Computer Science, pages 117–133. Springer, 2010.

[8]
S. Axler.
Measure, Integration & Real Analysis.
Graduate Texts in Mathematics. Springer International Publishing, 2019.

[9]
Christel Baier, Luca de Alfaro, Vojtech Forejt, and Marta Kwiatkowska.
Model checking probabilistic systems.
In Handbook of Model Checking, pages 963–999. Springer, 2018.

[10]
Christel Baier, Boudewijn R. Haverkort, Holger Hermanns, and Joost-Pieter Katoen.
Performance evaluation and model checking join forces.
Commun. ACM, 53(9):76–85, 2010.

[11]
Christel Baier and Joost-Pieter Katoen.
Principles of model checking.
MIT Press, 2008.

[12]
Jialu Bao, Nitesh Trivedi, Drashti Pathak, Justin Hsu, and Subhajit Roy.
Data-driven invariant learning for probabilistic programs.
In CAV (1), volume 13371 of Lecture Notes in Computer Science, pages 33–54. Springer, 2022.

[13]
Kevin Batz, Mingshuai Chen, Sebastian Junges, Benjamin Lucien Kaminski, Joost-Pieter Katoen, and Christoph Matheja.
Probabilistic program verification via inductive synthesis of inductive invariants.
In TACAS (2), volume 13994 of Lecture Notes in Computer Science, pages 410–429. Springer, 2023.

[14]
D. P. Bertsekas and S. E. Shreve.
Stochastic optimal control: The discrete-time case.
Athena Scientific, 1996.

[15]
Aleksandar Chakarov and Sriram Sankaranarayanan.
Probabilistic program analysis with martingales.
In CAV, volume 8044 of Lecture Notes in Computer Science, pages 511–526. Springer, 2013.

[16]
Ya-Chien Chang, Nima Roohi, and Sicun Gao.
Neural Lyapunov control.
In NeurIPS, pages 3240–3249, 2019.

[17]
Krishnendu Chatterjee, Hongfei Fu, and Amir Kafshdar Goharshady.
Termination analysis of probabilistic programs through positivstellensatz's.
In CAV (1), volume 9779 of Lecture Notes in Computer Science, pages 3–22. Springer, 2016.

[18]
Krishnendu Chatterjee, Hongfei Fu, Petr Novotný, and Rouzbeh Hasheminezhad.
Algorithmic analysis of qualitative and quantitative termination problems for affine probabilistic programs.
ACM Trans. Program. Lang. Syst., 40(2):7:1–7:45, 2018.

[19]
Krishnendu Chatterjee, Amir Kafshdar Goharshady, Tobias Meggendorfer, and Đorđe Žikelić.
Sound and complete certificates for quantitative termination analysis of probabilistic programs.
In CAV (1), volume 13371 of Lecture Notes in Computer Science, pages 55–78. Springer, 2022.

[20]
Krishnendu Chatterjee, Thomas A. Henzinger, Mathias Lechner, and Đorđe Žikelić.
A learner-verifier framework for neural network controllers and certificates of stochastic systems.
In TACAS (1), volume 13993 of Lecture Notes in Computer Science, pages 3–25. Springer, 2023.

[21]
Krishnendu Chatterjee, Petr Novotný, and Đorđe Žikelić.
Stochastic invariants for probabilistic termination.
In POPL, pages 145–160. ACM, 2017.

[22]
Fredrik Dahlqvist and Alexandra Silva.
Semantics of probabilistic programming: A gentle introduction.
In Gilles Barthe, Joost-Pieter Katoen, and Alexandra Silva, editors, Foundations of Probabilistic Programming, pages 1–42. Cambridge University Press, 2020.

[23]
Leonardo Mendonça de Moura and Nikolaj S. Bjørner.
Z3: an efficient SMT solver.
In TACAS, volume 4963 of Lecture Notes in Computer Science, pages 337–340. Springer, 2008.

[24]
Vojtech Forejt, Marta Z. Kwiatkowska, Gethin Norman, David Parker, and Hongyang Qu.
Quantitative multi-objective verification for probabilistic systems.
In TACAS, volume 6605 of Lecture Notes in Computer Science, pages 112–127. Springer, 2011.

[25]
Hongfei Fu and Krishnendu Chatterjee.
Termination of nondeterministic probabilistic programs.
In VMCAI, volume 11388 of LNCS, pages 468–490. Springer, 2019.

[26]
Timon Gehr, Sasa Misailovic, and Martin T. Vechev.
PSI: exact symbolic inference for probabilistic programs.
In CAV (1), volume 9779 of Lecture Notes in Computer Science, pages 62–83. Springer, 2016.

[27]
Timon Gehr, Samuel Steffen, and Martin T. Vechev.
\(\lambda\)PSI: exact inference for higher-order probabilistic programs.
In PLDI, pages 883–897. ACM, 2020.

[28]
Mirco Giacobbe, Daniel Kroening, and Julian Parsert.
Neural termination analysis.
In ESEC/SIGSOFT FSE, pages 633–645. ACM, 2022.

[29]
Andrew D. Gordon, Thomas A. Henzinger, Aditya V. Nori, and Sriram K. Rajamani.
Probabilistic programming.
In FOSE, pages 167–181. ACM, 2014.

[30]
O. Hernández-Lerma and J. B. Lasserre.
Discrete-time Markov control processes, volume 30 of Applications of Mathematics.
Springer-Verlag, 1996.

[31]
Mingzhang Huang, Hongfei Fu, Krishnendu Chatterjee, and Amir Kafshdar Goharshady.
Modular verification for almost-sure termination of probabilistic programs.
Proc. ACM Program. Lang., 3(OOPSLA):129:1–129:29, 2019.

[32]
O. Kallenberg.
Foundations of modern probability.
Springer Science & Business Media, 2006.

[33]
Dexter Kozen.
Semantics of probabilistic programs.
J. Comput. Syst. Sci., 22(3):328–350, 1981.

[34]
Marta Z. Kwiatkowska.
Quantitative verification: models techniques and tools.
In ESEC/SIGSOFT FSE, pages 449–458. ACM, 2007.

[35]
Marta Z. Kwiatkowska, Gethin Norman, and David Parker.
Quantitative analysis with the probabilistic model checker PRISM.
In QAPL, volume 153 of Electronic Notes in Theoretical Computer Science, pages 5–31. Elsevier, 2005.

[36]
Mathias Lechner, Đorđe Žikelić, Krishnendu Chatterjee, and Thomas A. Henzinger.
Stability verification in stochastic control systems via neural network supermartingales.
In AAAI, pages 7326–7336. AAAI Press, 2022.

[37]
Frederik Baymler Mathiesen, Simeon Craig Calvert, and Luca Laurenti.
Safety certification for stochastic systems via neural barrier functions.
IEEE Control. Syst. Lett., 7:973–978, 2023.

[38]
Annabelle McIver and Carroll Morgan.
Games, probability and the quantitative \(\mathrm{\mu}\)-calculus qm\(\mathrm{\mu}\).
In LPAR, volume 2514 of Lecture Notes in Computer Science, pages 292–310. Springer, 2002.

[39]
Annabelle McIver and Carroll Morgan.
Abstraction, Refinement and Proof for Probabilistic Systems.
Monographs in Computer Science. Springer, 2005.

[40]
Sean Meyn, Richard L. Tweedie, and Peter W. Glynn.
Markov Chains and Stochastic Stability.
Cambridge Mathematical Library. Cambridge University Press, 2 edition, 2009.

[41]
Carroll Morgan, Annabelle McIver, and Karen Seidel.
Probabilistic predicate transformers.
ACM Trans. Program. Lang. Syst., 18(3):325–353, 1996.

[42]
Van Chan Ngo, Quentin Carbonneaux, and Jan Hoffmann.
Bounded expectations: resource analysis for probabilistic programs.
In PLDI, pages 496–512. ACM, 2018.

[43]
Andrea Peruffo, Daniele Ahmed, and Alessandro Abate.
Automated and formal synthesis of neural barrier certificates for dynamical models.
In TACAS (1), volume 12651 of Lecture Notes in Computer Science, pages 370–388. Springer, 2021.

[44]
David Pollard.
A User's Guide to Measure Theoretic Probability.
Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2001.

[45]
H.L. Royden and P. Fitzpatrick.
Real Analysis.
Prentice Hall, 4 edition, 2010.

[46]
Armando Solar-Lezama.
Program Synthesis by Sketching.
PhD thesis, University of California at Berkeley, USA, 2008.

[47]
Armando Solar-Lezama, Liviu Tancau, Rastislav Bodík, Sanjit A. Seshia, and Vijay A. Saraswat.
Combinatorial sketching for finite programs.
In ASPLOS, pages 404–415. ACM, 2006.

[48]
Toru Takisaka, Yuichiro Oyabu, Natsuki Urabe, and Ichiro Hasuo.
Ranking and repulsing supermartingales for reachability in randomized programs.
ACM Trans. Program. Lang. Syst., 43(2):5:1–5:46, 2021.

[49]
Ilya Tkachev, Alexandru Mereacre, Joost-Pieter Katoen, and Alessandro Abate.
Quantitative model-checking of controlled discrete-time markov processes.
Inf. Comput., 253:1–35, 2017.

[50]
Franck van Breugel and James Worrell.
Towards quantitative verification of probabilistic transition systems.
In ICALP, volume 2076 of Lecture Notes in Computer Science, pages 421–432. Springer, 2001.

[51]
Di Wang, Jan Hoffmann, and Thomas W. Reps.
Central moment analysis for cost accumulators in probabilistic programs.
In PLDI, pages 559–573. ACM, 2021.

[52]
Peixin Wang, Hongfei Fu, Amir Kafshdar Goharshady, Krishnendu Chatterjee, Xudong Qin, and Wenjun Shi.
Cost analysis of nondeterministic probabilistic programs.
In PLDI, pages 204–220. ACM, 2019.