Technical University of Munich,
Germanykush.grover@in.tum.dehttps://orcid.org/0000-0003-4575-1302 Technical
University of Munich,
Germanyjan.kretinsky@in.tum.dehttps://orcid.org/0000-0002-8122-2881 Institute
of Science and Technology Austria, Vienna,
Austriatobias.meggendorfer@ist.ac.athttps://orcid.org/0000-0002-1712-2165
Technical University of Munich,
Germanymaxi.weininger@tum.dehttps://orcid.org/0000-0002-0163-2152 K. Grover,
J. Křetínský, T. Meggendorfer, and M. Weininger [500]Mathematics of computing
Markov processes [300]Mathematics of computing Continuous mathematics
[300]Computing methodologies Continuous models Bartek Klin, Sławomir Lasota,
and Anca Muscholl 3 33rd International Conference on Concurrency Theory
(CONCUR 2022) CONCUR 2022 CONCUR 2022 September 12–16, 2022 Warsaw, Poland 243
17

# Anytime Guarantees for Reachability in Uncountable Markov Decision Processes

Kush Grover  Jan Křetínský  Tobias Meggendorfer  Maximilian Weininger

###### Abstract

We consider the problem of approximating the reachability probabilities in
Markov decision processes (MDP) with uncountable (continuous) state and action
spaces. While there are algorithms that, for special classes of such MDP,
provide a sequence of approximations converging to the true value in the
limit, our aim is to obtain an algorithm with guarantees on the precision of
the approximation.

As this problem is undecidable in general, assumptions on the MDP are
necessary. Our main contribution is to identify sufficient assumptions that
are as weak as possible, thus approaching the ?boundary? of which systems can
be correctly and reliably analyzed. To this end, we also argue why each of our
assumptions is necessary for algorithms based on processing finitely many
observations.

We present two solution variants. The first one provides converging lower
bounds under weaker assumptions than typical ones from previous works
concerned with guarantees. The second one then utilizes stronger assumptions
to additionally provide converging upper bounds. Altogether, we obtain an
_anytime_ algorithm, i.e. yielding a sequence of approximants with known and
iteratively improving precision, converging to the true value in the limit.
Besides, due to the generality of our assumptions, our algorithms are very
general templates, readily allowing for various heuristics from literature in
contrast to, e.g., a specific discretization algorithm. Our theoretical
contribution thus paves the way for future practical improvements without
sacrificing correctness guarantees.

###### keywords:

Uncountable system, Markov decision process, discrete-time Markov control
process, probabilistic verification, anytime guarantee

##  1 Introduction

The standard formalism for modelling systems with both non-deterministic and
probabilistic behaviour are Markov decision processes (MDP) [47]. In the
context of many applications such as cyber-physical systems, states and
actions are used to model real-valued phenomena like position or throttle.
Consequently, the state space and the action space may be uncountably
infinite. For example, the intervals $[a,b]\times[c,d]\subseteq\mathbb{R}^{2}$
can model a safe area for a robot to move in or a set of available control
inputs such as acceleration and steering angle. This gives rise to MDP with
uncountable state- and action-spaces (sometimes called controlled discrete-
time Markov process [55, 56] or discrete-time Markov control process [14,
31]), with applications ranging from modelling a Mars rover [13, 26], over
water reservoir control [40] and warehouse storage management [42], to energy
control [55], and many more [45].

Although systems modelled by MDP are often safety-critical, the analysis of
uncountable systems is so complex that practical approaches for verification
and controller synthesis are usually based on ?best effort? learning
techniques, for example _reinforcement learning_. While efficient in practice,
these methods guarantee, even in the best case, convergence to the true result
only in the limit, e.g. [44], or for increasingly precise discretization, e.g.
[55, 35]. In line with the tradition of learning and to make the analysis more
feasible, the typical objectives considered for MDP are either finite-horizon
[41, 3] or discounted properties [21, 57, 28], together with restrictive
assumptions. Note that when it comes to approximation, discounted properties
effectively are finite-horizon. In contrast, ensuring safety of a reactive
system or a certain probability to satisfy its mission goals requires an
_unbounded_ horizon and reduces to optimizing the reachability probabilities.
Moreover, the safety-critical context requires _reliable_ bounds on the
probability, not an approximation with _unknown_ precision.

In this paper, we provide the first provably correct anytime algorithm for
(unbounded) reachability in uncountable MDP. As an _anytime_ algorithm, it can
at every step of the execution return correct lower and upper bounds on the
true value. Moreover, these bounds gradually converge to the true value,
allowing approximation up to an arbitrary precision. Since the problem is
undecidable, the core of our contribution is identifying sufficient conditions
on the uncountable MDP to allow for approximation.

Our _primary goal_ is to provide _conditions as weak as possible_ , thereby
pushing towards the boundary of which systems can be analyzed provably
correctly. To this end, we do not rely on any particular representation of the
system. Nonetheless, for classical scenarios, and, in particular, for finite
MDP, our conditions are mostly satisfied trivially.

Our _secondary goal_ is to derive the respective algorithms as an _extension
of value iteration_ (VI) [32, 47], while _avoiding drawbacks of
discretization_ -based approaches. VI is a de facto standard method for
numerical analysis of finite MDP, in particular with reachability objectives,
regarded as practically efficient and allowing for heuristics avoiding the
exploration of the complete state space, e.g. [12]. Interestingly, even for
finite MDP, anytime VI algorithms with precision guarantees are quite recent
[12, 22, 5, 48, 25]. Previous to that, the most used model checkers could
return arbitrarily wrong results [22]. Providing VI with precision guarantees
for general uncountable MDP is thus worthwhile on its own. Finally, while
discretization is conceptually simple, we prefer to provide a solution that
avoids the need to introduce arbitrary boundaries through gridding the whole
state space and, moreover, instead utilizes information from one ?cell? of the
grid in other places, too.

To summarize, while algorithmic aspects form an important motivation, our
primary contribution is theoretical: an explicit and complete set of generic
assumptions allowing for guarantees, disregarding practical efficiency at this
point. Consequently, while our approach lays foundations for further, more
tailored approaches, it is not to be seen as a competitor to the existing
practical, best-effort techniques, as these aim for a completely different
goal.

Our Contribution In this work, we provide the following:

Section 3:

A set of assumptions that allow for computing converging _lower_ bounds on the
reachability probability in MDP with uncountable state and action spaces. We
discuss in detail why they are weaker than usual, necessary, and applicable to
typically considered systems. With these assumptions, we extend the standard
(convergent but precision-ignorant) VI to this general setting.

Section 4:

An additional set of assumptions that yield the first _anytime_ algorithm,
i.e. with provable bounds on the precision/error of the result, converging to
0. We combine the preceding algorithm with the technique of _bounded real-time
dynamic programming (BRTDP)_ [43] and provide also converging _upper_ bounds
on the reachability probability.

Section 5:

A discussion of theoretical extensions and practical applications.

Related work For detailed _theoretical_ treatment of reachability and related
problems on uncountable MDP, see e.g. [56, 14]. Reachability on uncountable
MDP generalizes numerous problems known to be undecidable. For example, we can
encode the halting problem of (probabilistic) Turing machines by encoding the
tape content as real value. Similarly, almost-sure termination of
probabilistic programs (undecidable [36]) is a special case of reachability on
general uncountable MDP (see e.g. [19]). As precise reachability analysis is
undecidable even for non-stochastic linear hybrid systems [29], many works
turn their attention to more relaxed notions such as $\delta$-reachability,
e.g. [52], and/or employ many assumptions.

In order to obtain precision bounds, we assume that the _value_ function,
mapping states to their reachability probability, is _Lipschitz continuous_
(and that we know the Lipschitz constant). This is slightly _weaker_ than the
classical approach of assuming Lipschitz continuity of the _transition_
function (and knowledge of the constant), e.g. [2, 53]. In particular, these
assumptions (i) imply our assumption (as we show in Section B.2.1) and (ii)
are used even in the simpler settings of finite-horizon and discounted reward
scenarios [7, 2, 53, 55] or even more restricted settings to obtain practical
efficiency, e.g. [39]. In contrast to our approach, they are not anytime
algorithms and require treatment of the whole state space.

To provide context, we outline how continuity is used (explicitly or
implicitly) in related work and mention their respective results. Firstly,
[28, 51] assume Lipschitz continuity, _but not_ explicit knowledge of the
constant. In essence, these approaches solve the problem by successively
increasing internal parameters.The parameters then eventually cross a bound
implied by the Lipschitz constant, yielding an ?eventual correctness?. In
particular, they provide ?convergence in the limit? or ?probably approximately
correct? results, but no bounds on the error or the convergence rate; these
would depend on knowledge of the constant.

Secondly, [21, 44, 2, 53, 55] (and our work) assume Lipschitz continuity _and_
knowledge of the constant. Relying on the constant being provided externally,
these works derive guarantees. Previously, the guarantees given are weaker
than our convergent anytime bounds: Either convergence in the limit [44] or a
bound on a discretization error, relativized to sub-optimal strategies [21] or
bounded horizon [2, 53, 55].

Several of the above mentioned works employ _discretization_ [21, 2, 53, 55].
This method is quite general, but obtaining any bounds on the error requires
continuity assumptions [1]. Further, there are works that use other
assumptions: [27, 26] use _reinforcement learning_ methods to tackle
reachability and more general problems, without any continuity assumption.
However, they do not provide any guarantees. See [57] for a detailed
exposition of similar approaches. Assuming an abstraction is given,
_abstraction and bisimulation_ approaches, e.g. [24, 23], provide guarantees,
but only on the lower bounds. With significant assumptions on the system’s
structure, _symbolic_ approaches [41, 59, 49, 17] may even obtain exact
solutions.

##  2 Preliminaries

In this section, we recall basics of probabilistic systems and set up the
notation. As usual, $\mathbb{N}$ and $\mathbb{R}$ refer to the (positive)
natural numbers and real numbers, respectively. For a set $S$,
$\mathbbm{1}_{S}$ denotes its _characteristic function_ , i.e.
$\mathbbm{1}_{S}(x)=1$ if $x\in S$ and $0$ otherwise. We write $S^{\star}$ and
$S^{\omega}$ to refer to the set of finite and infinite sequences comprising
elements of $S$, respectively.

We assume familiarity with basic notions of measure theory, e.g.  _measurable
set_ or _measurable function_ , as well as probability theory, e.g.
_probability spaces_ and _measures_ [11]. For a measure space $X$ with sigma-
algebra $\Sigma_{X}$, $\Pi(X)$ denotes the set of all probability measures on
$X$. For a measure $\mu\in\Pi(X)$, we write $\mu(Y)=\int\mathbbm{1}_{Y}\,d\mu$
to denote the mass of a measurable set $Y\in\Sigma_{X}$ (also called _event_).
For two probability measures $\mu$ and $\nu$, the _total variation distance_
is defined as
$\delta_{TV}(\mu,\nu):=2\cdot\sup_{Y\in\Sigma_{X}}\lvert\mu(Y)-\nu(Y)\rvert$.
Some event happens _almost surely_ (a.s.) w.r.t. some measure $\mu$ if it
happens with probability $1$. We write $\operatorname{supp}(\mu)$ to denote
the _support_ of the probability measure $\mu$.

######  Remark 2.1.

It is surprisingly difficult to give a well-defined notion of support for
measures in general. Intuitively, $\operatorname{supp}(\mu)$ describes the
?smallest? set which $\mu$ assigns a value of $1$. However, this is not well-
defined for general measures. We discuss these issues and a proper definition
in Appendix E. Throughout this work, similar subtle issues related to measure
theory arise. For the sake of readability, these are mostly delegated to
footnotes or the appendix, and readers may safely skip over these points.

We work with _Markov decision processes_ (MDP) [47], a widely used model to
capture both non-determinism and probability. We consider uncountable state
and action spaces.

######  Definition 2.2.

A _(continuous-space, discrete-time) Markov decision process (MDP)_ is a tuple
$\mathcal{M}=(S,Act,{Av},\Delta)$, where $S$ is a compact set of _states_
(with topology $\mathcal{T}_{S}$ and Borel $\sigma$-algebra
$\Sigma_{S}=\mathfrak{B}(\mathcal{T}_{S})$), $Act$ is a compact set of
_actions_ (with topology $\mathcal{T}_{Act}$ and Borel $\sigma$-algebra
$\Sigma_{Act}=\mathfrak{B}(\mathcal{T}_{Act})$), ${Av}\colon
S\to\Sigma_{Act}\setminus\\{\emptyset\\}$ assigns to every state a non-empty,
measurable, and compact set of _available actions_ , and $\Delta\colon S\times
Act\to\Pi(S)$ is a _transition function_ that for each state $s$ and
(available) action $a\in{Av}(s)$ yields a probability measure over successor
states (i.e. a Markov Kernel). An MDP is called _finite_ if $\lvert
S\rvert<\infty$ and $\lvert Act\rvert<\infty$.

See [47, Sec. 2.3] and [8, Chp. 9] for a more detailed discussion on the
technical considerations arising from uncountable state and action spaces.
Note that we assume the set of available actions to be non-empty. This means
that the system can never get ?stuck? in a degenerate state without
successors. _Markov chains_ are a special case of MDP where
$\lvert{Av}(s)\rvert=1$ for all $s\in S$, i.e. a completely probabilistic
system without any non-determinism. Our presented methods thus are directly
applicable to Markov chains as well.

Given a measure $\mu\in\Pi(X)$ and a measurable function $f\colon
X\to\mathbb{R}$ mapping elements of a set $X$ to real numbers, we write
$\mu\langle f\rangle:=\int f(x)\,d\mu(x)$ to denote the integral of $f$ with
respect to $\mu$. For example, $\Delta(s,a)\langle f\rangle$ denotes the
expected value $\mathbb{E}_{s^{\prime}\sim\Delta(s,a)}f(s^{\prime})$ of
$f\colon S\to\mathbb{R}$ over the successors of $s$ under action $a$.
Moreover, abusing notation, for some set of state $S^{\prime}\subseteq S$ and
function ${Av}^{\prime}\colon S^{\prime}\to Act$, we write
$S^{\prime}\times{Av}^{\prime}=\\{(s,a)\mid s\in
S^{\prime},a\in{Av}^{\prime}(s)\\}$ to denote the set of state-action pairs
with states from $S^{\prime}$ under ${Av}^{\prime}$.

An _infinite path_ in an MDP is some infinite sequence
$\rho=s_{1}a_{1}s_{2}a_{2}\cdots\in(S\times{Av})^{\omega}$, such that for
every $i\in\mathbb{N}$ we have
$s_{i+1}\in\operatorname{supp}(\Delta(s_{i},a_{i}))$. A _finite path_ (or
_history_) $\varrho=s_{1}a_{1}s_{2}a_{2}\dots
s_{n}\in(S\times{Av})^{\star}\times S$ is a non-empty, finite prefix of an
infinite path of length $\lvert\varrho\rvert=n$, ending in state $s_{n}$,
denoted by $last(\varrho)$. We use $\rho(i)$ and $\varrho(i)$ to refer to the
$i$-th state in an (in)finite path. We refer to the set of finite (infinite)
paths of an MDP $\mathcal{M}$ by $\mathsf{FPaths}_{\mathcal{M}}$
($\mathsf{Paths}_{\mathcal{M}}$). Analogously, we write
$\mathsf{FPaths}_{\mathcal{M},s}$ ($\mathsf{Paths}_{\mathcal{M},s}$) for all
(in)finite paths starting in $s$.

In order to obtain a probability measure, we first need to eliminate the non-
determinism. This is done by a so-called _strategy_ (also called _policy_ ,
_controller_ , or _scheduler_). A strategy on an MDP
$\mathcal{M}=(S,Act,{Av},\Delta)$ is a function
$\pi\colon\mathsf{FPaths}_{\mathcal{M}}\to\Pi(Act)$, s.t.
$\operatorname{supp}(\pi(\varrho))\subseteq{Av}(last(\varrho))$. The set of
all strategies is denoted by $\Pi_{\mathcal{M}}$. Intuitively, a strategy is a
?recipe? describing which step to take in the current state, given the
evolution of the system so far.

Given an MDP $\mathcal{M}$, a strategy $\pi\in\Pi_{\mathcal{M}}$, and an
initial state $s_{0}$, we obtain a measure on the set of infinite paths
$\mathsf{Paths}_{\mathcal{M}}$, which we denote as
$\mathsf{Pr}_{\mathcal{M},s_{0}}^{\pi}$. See [47, Sec. 2] for further details.
Thus, given a measurable set $A\subseteq\mathsf{Paths}_{\mathcal{M}}$, we can
define its maximal probability starting from state $s_{0}$ under any strategy
by
$\mathsf{Pr}_{\mathcal{M},s_{0}}^{\sup}[A]:={\sup}_{\pi\in\Pi_{\mathcal{M}}}\mathsf{Pr}_{\mathcal{M},s_{0}}^{\pi}[A].$
Depending on the structure of $A$ it may be the case that no optimal strategy
exists and we have to resort to the supremum instead of the maximum. This may
already arise for finite MDP, see [15].

For an MDP $\mathcal{M}=(S,Act,{Av},\Delta)$ and a set of _target states_
$T\subseteq S$, _(unbounded) reachability_ refers to the set $\Diamond
T=\\{\rho\in\mathsf{Paths}_{\mathcal{M}}\mid\exists
i\in\mathbb{N}.\leavevmode\nobreak\ \rho(i)\in T\\}$, i.e. all paths which
eventually reach $T$. The set $\Diamond T$ is measurable if $T$ is measurable
[55, Sec. 3.1], [56, Sec. 2].

Now, it is straightforward to define the _maximal reachability problem_ of a
given set of states. Given an MDP $\mathcal{M}$, target set $T$, and state
$s_{0}$, we are interested in computing the maximal probability of eventually
reaching $T$, starting in state $s_{0}$. Formally, we want to compute the
_value_ of the state $s_{0}$, defined as
$\mathcal{V}(s_{0}):=\mathsf{Pr}_{\mathcal{M},s_{0}}^{\sup}[\Diamond
T]={\sup}_{\pi\in\Pi_{\mathcal{M}}}\mathsf{Pr}_{\mathcal{M},s_{0}}^{\pi}[\Diamond
T].$ This state value function satisfies a straightforward fixed point
equation, namely

$\mathcal{V}(s)=1\quad\text{if $s\in
T$}\qquad\mathcal{V}(s)={\sup}_{a\in{Av}(s)}\Delta(s,a)\langle\mathcal{V}\rangle\qquad\text{otherwise.}$
(1)

Moreover, $\mathcal{V}$ is the _smallest_ fixed point of this equation [8,
Prop. 9.8, 9.10], [56, Thm. 3]. In our approach, we also deal with values of
state-action pairs $(s,a)\in S\times{Av}$, where
$\mathcal{V}(s,a):=\Delta(s,a)\langle\mathcal{V}\rangle$. Intuitively, this
represents the value achieved by choosing action $a$ in state $s$ and then
moving optimally. Clearly, we have that
$\mathcal{V}(s)=\sup_{a\in{Av}(s)}\mathcal{V}(s,a)$. See [18, Sec. 4] for a
discussion of reachability on finite MDP and [56] for the general case.

In this work, we are interested in _approximate solutions_ due to the
following two reasons. Firstly, obtaining precise solutions for MDP is
difficult already under strict assumptions and undecidable in our general
setting.(1)(1)(1)For example, one can encode the tape of a Turing machine into
the binary representation of a real number and reduce the halting problem to a
reachability query. We thus resort to approximation, allowing for much lighter
assumptions. Secondly, by considering approximation we are able to apply many
different optimization techniques, potentially leading to algorithms which are
able to handle real-world systems, which are out of reach for precise
algorithms even for finite MDP [12].

We are interested in two types of approximations. Firstly, we consider
approximating the value function in the limit, without knowledge about how
close we are to the true value. This is captured by a semi-decision procedure
for queries of the form $\mathsf{Pr}_{\mathcal{M},s}^{\sup}[\Diamond T]>\xi$
for a threshold $\xi\in[0,1]$. We call this problem ApproxLower. Secondly, we
consider the variant where we are given a precision requirement
$\varepsilon>0$ and obtain $\varepsilon$-optimal values $(l,u)$, i.e. values
with $\mathcal{V}(s_{0})\in[l,u]$ and $0\leq u-l<\varepsilon$. We refer to
this variant as ApproxBounds.

##  3 Converging Lower Bounds

In this section, we present the first set of assumptions, enabling us to
compute _converging lower bounds_ on the true value, solving the ApproxLower
problem. In Section 3.1, we discuss each assumption in detail and argue on an
intuitive level why it is necessary by means of counterexamples. With the
assumptions in place, in Section 3.2 we then present our first algorithm, also
introducing several ideas we employ again in the following section.

Our assumptions and algorithms are motivated by _value iteration_ (VI) [32],
which we briefly outline. In a nutshell, VI boils down to repeatedly applying
an iteration operator to a value vector $v_{n}$. For example, the canonical
value iteration for reachability on finite MDP starts with $v_{0}(s)=1$ for
all $s\in T$ and $0$ otherwise and then iterates

$v_{n+1}(s)={\max}_{a\in Act(s)}{\sum}_{s^{\prime}\in
S}\Delta(s,a,s^{\prime})\cdot v_{n}(s^{\prime})$ (2)

for all $s\notin T$. The vector $v_{n}$ converges monotonically from below to
the true value for all states. We mention two important points. Firstly, the
iteration can be applied ?asynchronously?. Instead of updating all states in
every iteration, we can pick a single state and only update its value. The
values $v_{n}$ still converge to the correct value as long as all states are
updated infinitely often. Secondly, instead of storing a value per state, we
can store a value for each state-action pair and obtain the state value as the
maximum of these values. Both points are a technical detail for finite MDP,
however they play an essential role in our uncountable variant. See Section
A.1 for more details on VI for finite MDP.

In the uncountable variant of Equation 2, $v$ is a function, $Act(s)$ is
potentially uncountable, and the sum is replaced by integration. As in this
setting the problem is undecidable, naturally we have to employ some
assumptions. Our goal is to sufficiently imitate the essence of Equation 2,
obtaining convergence without being overly restrictive. In particular, we want
to (i) represent (an approximation of) $v_{n}$ using finite memory, (ii)
safely approximate the maximum and integration, and (iii) select appropriate
points to update $v_{n}$.

###  3.1 Assumptions

Before discussing each assumption in detail, we first put them into context.
As we argue in the following, most of our assumptions typically hold
implicitly. Still, by stating even basic computability assumptions in a form
as weak as possible, we avoid ?hidden? assumptions, e.g. by assuming that the
state space is a subset of $\mathbb{R}^{d}$. Two of our assumptions are more
restrictive, namely Assumption C: Value Lipschitz Continuity (Section 3.1.3)
and, introduced later, Assumption D: Absorption (Section 4.1.2). However, they
are also often used in related works, as we detail in the respective sections.
Moreover, in light of previous results, the necessity of restrictive
assumptions is to be expected: Computing bounds is hard or even undecidable
already for very restricted classes. Aside from the discussion in the
introduction, we additionally mention two further cases. In the setting of
probabilistic programs (which are a very special case of uncountable MDP),
deciding almost sure termination for a fixed initial state (which is a
severely restricted subclass of reachability on uncountable MDP without non-
determinism) is an actively researched topic with recent advances, see e.g.
[33, 34], and shown to be $\Pi_{2}^{0}$-complete [36], i.e. highly
undecidable. In [30] and the references therein, the authors present
(un-)decidability results for hybrid automata, which are a special case of
uncountable MDP _without any stochastic dynamics_ (flow transitions can be
modelled as actions indicating the delay). As such, it is to be expected that
the general class of models we consider has to be pruned very strictly in
order to hope for any decidability results.

######  Remark 3.1.

As already mentioned, we want to provide assumptions which are as general as
possible. Importantly, we avoid (unnecessarily) assuming any particular
representation of the system. Our motivation is to ultimately identify the
boundary of what is necessary to derive guarantees. While our assumptions are
motivated by VI and built around Equation 2, we note that being able to
represent the state values and evaluate (some aspect of) the transition
dynamics intuitively are a necessity for _any_ method dealing with such
systems. We do not claim that our framework of assumptions is the only way to
approach the problem, instead we provide arguments why it is a sensible way to
do so.

####  3.1.1 A: Basic Assumptions (Asm. A1-A4)

We first present a set of basic computability assumptions (A1-A4). These are
essential, since for uncountable systems even the simplest computations are
intractable without any assumptions. More specifically, such systems cannot be
given explicitly (due to their infinite size), but instead have to be
described symbolically by, e.g., differential equations. Thus, we necessarily
require some notion of computability and structural properties for each part
of this symbolic description. And indeed, each assumption essentially
corresponds to one part of the MDP description (Metric Space to $S\times Act$,
Maximum Approximation to ${Av}$, Transition Approximation to $\Delta$, and
Target Computability to $T$). They are weak and hold on practically all
commonly considered systems (see Section B.1). In particular, finite MDP and
discrete components are trivially subsumed by considering the discrete metric.

A1: Metric Space

$S$ and $Act$ are metric spaces with (computable) metrics
$\operatorname{d}_{S}$ and $\operatorname{d}_{Act}$, respectively, and
$\operatorname{d}_{\times}$ is a compatible(2)(2)(2)For two pairs $(s,a)$ and
$(s,a^{\prime})$ we have that
$k\cdot\operatorname{d}_{Act}(a,a^{\prime})\leq\operatorname{d}_{\times}((s,a),(s,a^{\prime}))\leq
K\cdot\operatorname{d}_{Act}(a,a^{\prime})$ for some constants $k,K\geq 0$,
analogous for $\operatorname{d}_{S}$, achieved by, e.g.
$\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))\coloneqq\operatorname{d}_{S}(s,s^{\prime})+\operatorname{d}_{Act}(a,a^{\prime})$.
metric on the space of state-action pairs $S\times{Av}$,

A2: Maximum Approximation

For each state $s$ and computable Lipschitz $f:{Av}(s)\to[0,1]$, the value
$\max_{a\in{Av}(s)}f(a)$ can be under-approximated to arbitrary precision.

A3: Transition Approximation

For each state-action pair $(s,a)$ and Lipschitz $g:S\to[0,1]$ which can be
under-approximated to arbitrary precision, the successor expectation
$\Delta(s,a)\langle g\rangle$ can be under-approximated to arbitrary
precision.

A4: Target Computability

The target set $T$ is decidable, i.e. we are given a computable predicate
which, given a state $s$, decides whether $s\in T$.

We denote the approximations for A2 and A3 by $\textsc{Approx}_{\leq}$, i.e.
given a pair $(s,a)$ and functions $f$, $g$ as in the assumptions, we write
(abusing notation)
$\textsc{Approx}_{\leq}(\max_{a\in{Av}(s)}f(a),\varepsilon)$ and
$\textsc{Approx}_{\leq}(\Delta(s,a)\langle g\rangle,\varepsilon)$ for
approximation of the respective values up to precision $\varepsilon$, i.e.
$0\leq\max_{a\in{Av}(s)}f(a)-\textsc{Approx}_{\leq}(\max_{a\in{Av}(s)}f(a),\varepsilon)\leq\varepsilon$
and analogous for $\Delta(s,a)\langle g\rangle$. Note that A2 and A3 are
satisfied if we can sample densely in ${Av}(s)$ and approximate $\Delta(s,a)$.

####  3.1.2 B: Sampling (Asm. B.VI)

As there are uncountably many states, we are unable to explicitly update all
of them at once and instead update values asynchronously. Moreover, as there
may also be uncountably many actions, we instead store and update the values
of state-action pairs. Together, we need to pick state-action pairs to update.
We delegate this choice to a selection mechanism GetPair, an oracle for state-
action pairs. We allow for GetPair to be ?stateful?, i.e. the sampled state-
action pair may depend on previously returned pairs. This is required in, for
example, round-robin or simulation-based approaches. We only require a basic
notion of fairness in order to guarantee that we do not miss out on any
information. Note the additional identifier .VI (_value iteration_) on the
assumption name; later on, a similar, but weaker variant (B.BRTDP) is
introduced.

B.VI: State-Action Sampling

Let
$S^{\Diamond}=\\{last(\varrho)\mid\varrho\in\mathsf{FPaths}_{\mathcal{M},s}\\}$
the set of all reachable states. Then, for any $\varepsilon>0$, $s\in
S^{\Diamond}$, and $a\in{Av}(s)$ we have that GetPair eventually yields a pair
$(s^{\prime},a^{\prime})$ with
$\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))<\varepsilon$ and
$\delta_{TV}(\Delta(s,a),\Delta(s^{\prime},a^{\prime}))<\varepsilon$
a.s.(3)(3)(3)Technically, it is sufficient to satisfy this property on any
subset of $S^{\Diamond}$ which only differs from it up to measure 0. More
precisely, we only require that this assumption holds for
$S^{\Diamond}=\operatorname{supp}(\mathsf{Pr}_{\mathcal{M},s}^{\sup})$, i.e.
the set of all reachable paths with non-zero measure. We omit this rather
technical notion and the discussion it entails in order to avoid distracting
from the central results of this work.

Essentially, this means that GetPair provides a way to ?exhaustively? generate
all behaviours of the system up to a precision of $\varepsilon$. This fairness
assumption is easily satisfied under usual conditions. For example, if
$S\times{Av}$ is a bounded subset of $\mathbb{R}^{d}$, we can randomly sample
points in that space or consider increasingly dense grids. Alternatively, if
we can sample from the set of actions and from the distributions of $\Delta$,
GetPair can be implemented by sampling paths of random length, following
random actions. Note that we can view the procedure as a ?template?: Instead
of requiring a concrete method to acquire pairs to update, we leave this open
for generality; we discuss implications of this in Sections 5.1 and 5.3.

The requirement on total variation may seem unnecessary, especially given that
we will also assume continuity. However, otherwise we could, for example, miss
out on solitary actions which are the ?witnesses? for a state’s value: suppose
that ${Av}(s)=[0,1]$ and $\Delta(s,0)$ moves to the goal, while $\Delta(s,a)$
just loops back to $s$. Only selecting actions close to $a=0$ w.r.t. the
product metric is not sufficient to observe that we can move to the goal. Note
that this would not be necessary if we assumed continuity of the transition
function – selecting ?nearby? actions then also yields ?similar? behaviour.

####  3.1.3 C: Lipschitz Continuity

Finally, we present our already advertised continuity assumption. For
simplicity, we give it in its strict form and discuss relaxations later in
Section 5.2. Intuitively, Lipschitz continuity allows us to extrapolate the
behaviour of the system from a single state to its surroundings.

C: Value Lipschitz Continuity

The value functions $\mathcal{V}(s)$ and $\mathcal{V}(s,a)$ are Lipschitz
continuous with _known_ constants $C_{S}$ and $C_{\times}$, i.e. for all
$s,s^{\prime}\in S$ and $a\in{Av}(s),a^{\prime}\in{Av}(s^{\prime})$ we have

$\displaystyle\lvert\mathcal{V}(s)-\mathcal{V}(s^{\prime})\rvert$
$\displaystyle\leq C_{S}\cdot\operatorname{d}_{S}(s,s^{\prime})$
$\displaystyle\lvert\mathcal{V}(s,a)-\mathcal{V}(s^{\prime},a^{\prime})\rvert$
$\displaystyle\leq
C_{\times}\cdot\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))$

This requirement may seem quite restrictive at first glance. Indeed, it is the
only one in this section to not usually hold on ?standard? systems. However,
in order to obtain any kind of (provably correct) bounds, some notion of
continuity is elementary, since otherwise we cannot safely extrapolate from
finitely many observations to an uncountable set. The immediately arising
questions are (i) why Lipschitz continuity is necessary compared to, e.g.,
regular or uniform continuity, and (ii) why knowledge of the Lipschitz
constant is required. For the first point, note that we want to be able to
extrapolate from values assigned to a single state to its immediate
surroundings. While continuity means that the values in the surroundings do
not ?jump?, it does not give us any way of bounding the rate of change, and
this rate may grow arbitrarily (for example, consider the continuous but not
Lipschitz function $\sin(\frac{1}{x})$ for $x>0$). So, also relating to the
second point, without knowledge of the Lipschitz constant, regular continuity
and Lipschitz continuity are (mostly) equivalent from a computational
perspective: The function does not have discontinuities, but we cannot safely
estimate the rate of change in general. To illustrate this point further, we
give an intuitive example.

######  Example 3.2.

We construct an MDP with a periodic, Lipschitz continuous value function, as
illustrated in Figure 1 and formally defined below. Intuitively, for a given
period width $w$ (e.g. 0.25) and a periodic function $f$ (e.g. a triangle
function), a state $s$ between $0$ and $w$ moves to a target or sink with
probability $f(s)$. All larger states $s\geq w$ transition to $s-w$ with
probability 1. The value function thus is periodic and Lipschitz continuous,
see Figure 1 for a possible value function and Section B.2.3 for a formal
definition.

0.250.50.75$0$$0.5$$1$$S$$\mathcal{V}(s)$ Figure 1:  The value function of
Example 3.2, showing that knowledge of the constant is important.

For a finite number of samples, we can choose $f$ and $w$ such that all
samples achieve a value of $1$. Nevertheless, we cannot conclude anything
about states we have not sampled yet: _Without knowledge of the constant, we
cannot extrapolate from samples_.

We note the underlying connection to the Nyquist-Shannon sampling theorem [50,
Thm. 1]. Intuitively, the theorem states that, for a function that contains no
frequencies higher than $W$, it is completely determined by giving its
ordinates at a series of points spaced $0.5\cdot W$ apart. If we know the
Lipschitz constant, this gives us a way of bounding the ?frequency? of the
value function, and thus allows us to determine it by sampling a finite number
of points. On the other hand, without the Lipschitz constant, we do not know
the frequency and cannot judge whether we are ?undersampling?.

Since we do not assume any particular representation of the transition system,
we cannot derive such constants in general. Instead, these would need to be
obtained by, e.g., domain knowledge, or tailored algorithms. As in previous
approaches [21, 44, 2, 53, 55], we thus resort to assuming that we are given
this constant, offloading this (highly non-trivial) step. Recall that
Lipschitz continuity of the transition function implies Lipschitz continuity
of the value function (see Section B.2.1), but can potentially be checked more
easily.

###  3.2 Assumptions Applied: Value Iteration Algorithm

Before we present our new algorithm, we explain how our assumptions allow us
to lift VI to the uncountable domain. Contrary to the finite state setting, we
are unable to store precise values for each state explicitly, since there are
uncountably many states. Hence, the algorithm exploits the Lipschitz-
continuity of the value function as follows. Assume that we know that the
value of a state $s$ is bounded from below by a value $l$, i.e.
$\mathcal{V}(s)\geq l$. Then, by Lipschitz-continuity of $\mathcal{V}$, we
know that the value of a state $s^{\prime}$ is bounded by
$l-\operatorname{d}_{S}(s,s^{\prime})\cdot C_{S}$. More generally, if we are
given a finite set of states $\mathsf{Sampled}$ with correct lower bounds
$\widehat{\mathsf{L}}\colon\mathsf{Sampled}\to[0,1]$, we can safely extend
these values to the whole state space by

$\mathsf{L}(s)\coloneqq{\max}_{s^{\prime}\in\mathsf{Sampled}}\left(\widehat{\mathsf{L}}(s^{\prime})-C_{S}\cdot\operatorname{d}_{S}(s,s^{\prime})\right).$

Since $\mathcal{V}(s)\geq\widehat{\mathsf{L}}(s)$ for all
$s\in\mathsf{Sampled}$, we have $\mathcal{V}(s)\geq\mathsf{L}(s)$ for all
$s\in S$, i.e. $\mathsf{L}(\cdot)$ is a valid lower bound. We thus obtain a
lower bound for all of the uncountably many states, described symbolically as
a combination of finitely many samples. See Figure 2 for an illustration.

0.511.5$0$$0.2$$0.4$$S$ Figure 2: Example of the function extension on the set
$[0,2]$ with a Lipschitz constant of $C_{S}=1$. Dots represent stored values
in $\widehat{\mathsf{L}}$, while the solid line represents the extrapolated
function $\mathsf{L}$. Note that it is possible to have
$\widehat{\mathsf{L}}(s)<\mathsf{L}(s)$, as seen in the graph.

This is sufficient to deal with Markov chains, but for MDPs we additionally
need to take care of the (potentially uncountably many) actions. Recall that
value iteration updates state values with the maximum over available actions,
$v_{n+1}(s)=\max_{a\in{Av}(s)}\Delta(s,a)\langle v_{n}\rangle$. This is
straightforward to compute when there are only finitely many actions, but in
the uncountable case obtaining
$\mathsf{L}(s)=\sup_{a\in{Av}(s)}\mathsf{L}(s,a)$ is much more involved. We
apply the idea of Lipschitz continuity again, storing values for a set
$\mathsf{Sampled}$ of state-action pairs instead of only states. We bound the
value of every state-action pair by

$\mathsf{L}(s,a)\coloneqq{\max}_{(s^{\prime},a^{\prime})\in\mathsf{Sampled}}\left(\widehat{\mathsf{L}}(s^{\prime},a^{\prime})-\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))\cdot
C_{\times}\right)$ (3)

Observe that $\mathsf{L}(s,a)$ is computable and Lipschitz-continuous as well,
so by Maximum Approximation we can approximate the bound of any state, i.e.
$\mathsf{L}(s)=\max_{a\in{Av}(s)}\mathsf{L}(s,a)$, based on such a finite set
of values assigned to state-action pairs. (Recall that ${Av}(s)$ is compact
and $\mathsf{L}(s,a)$ continuous, hence the maximum is attained.)
Consequently, we can also under-approximate
$\Delta(s,a)\langle\mathsf{L}\rangle$ by Transition Approximation. To avoid
clutter, we omit the following two special cases in the definition of
$\mathsf{L}(s,a)$: Firstly, if $\mathsf{Sampled}=\emptyset$, we naturally set
$\mathsf{L}(s,a)=0$. Secondly, if all pairs $(s^{\prime},a^{\prime})$ are too
far away for a sensible estimate, i.e. if Equation 3 was yielding
$\mathsf{L}(s,a)<0$, we also set $\mathsf{L}(s,a)$ to $0$.

Algorithm 1 The Value Iteration (VI) Algorithm for MDPs with general state-
and action-spaces.

1:ApproxLower query with threshold $\xi$, satisfying A1–A4, B.VI and C.

2:yes, if $\mathcal{V}(s_{0})>\xi$.

3:$\mathsf{Sampled}\leftarrow\emptyset,\mathsf{t}\leftarrow 1$
$\triangleright$ Initialize

4:while
$\textsc{Approx}_{\leq}(\mathsf{L}(s_{0}),\textsc{Precision}(\mathsf{t}))\leq\xi$
do

5: $(s,a)\leftarrow\textsc{GetPair}$ $\triangleright$ Sample state-action pair

6: if $s\in T$ then $\widehat{\mathsf{L}}(s,\cdot)\leftarrow 1$
$\triangleright$ Handle target states

7: else
$\widehat{\mathsf{L}}(s,a)\leftarrow\textsc{Approx}_{\leq}(\Delta(s,a)\langle\mathsf{L}\rangle,\textsc{Precision}(\mathsf{t}))$
$\triangleright$ Update $\widehat{\mathsf{L}}$

8: $\mathsf{Sampled}\leftarrow\mathsf{Sampled}\mathbin{\cup}\\{(s,a)\\}$,
$\mathsf{t}\leftarrow\mathsf{t}+1$

9:return yes

We present VI for MDPs with general state- and action-spaces in Algorithm 1.
It depends on $\textsc{Precision}(\mathsf{t})$, a sequence of precisions
converging to zero in the limit, e.g.
$\textsc{Precision}(\mathsf{t})=\frac{1}{\mathsf{t}}$. The algorithm executes
the main loop until the current approximation of the lower bound of the
initial state $\mathsf{L}(s_{0})=\max_{a\in{Av}(s_{0})}\mathsf{L}(s_{0},a)$
exceeds the given threshold $\xi$. Inside the loop, the algorithm updates
state-action pairs yielded by GetPair. For target states, the lower bound is
set to $1$. Otherwise, we set the bound of the selected pair to an
approximation of the expected value of $\mathsf{L}$ under the corresponding
transition. Here is the crucial difference to VI in the finite setting:
Instead of using Equation 2, we have to use Equation 3 and
$\textsc{Approx}_{\leq}$, the approximations that exist by assumption, see
Section 3.1.1. Since $\textsc{Precision}(\mathsf{t})$ converges to zero, the
approximations eventually get arbitrarily fine. The procedure
$\textsc{Precision}(\mathsf{t})$ may be adapted heuristically in order to
speed up computation. For example, it may be beneficial to only approximate up
to $0.01$ precision at first to quickly get a rough overview. We show that
Algorithm 1 is correct, i.e. the stored values (i) are lower bounds and (ii)
converge to the true values in Section E.1. Here, we only provide a sketch,
illustrating the main steps.

######  Theorem 3.3.

Algorithm 1 is correct under Assumptions A1–A4, B.VI, and C, i.e. it outputs
yes iff $\mathcal{V}(s)>\xi$.

######  Proof 3.4 (Proof sketch).

First, we show that
$\mathsf{L}_{\mathsf{t}}(s)\leq\mathsf{L}_{\mathsf{t}+1}(s)\leq\mathcal{V}(s)$
by simple induction on the step. Initially, we have $\mathsf{L}_{1}(s)=0$,
obviously satisfying the condition. The updates in Lines 6 and 7 both keep
correctness, i.e. $\mathsf{L}_{\mathsf{t}+1}(s)\leq\mathcal{V}(s)$, proving
the claim.

Since $\mathsf{L}_{\mathsf{t}}$ is monotone as argued above, its limit for
$\mathsf{t}\to\infty$ is well defined, denoted by $\mathsf{L}_{\infty}$. By
State-Action Sampling, the set of accumulation points of $s_{\mathsf{t}}$
contains all reachable states $S^{\Diamond}$. We then prove that
$\mathsf{L}_{\infty}$ satisfies the fixed point equation Equation 1. For this,
we use the second part of the assumption on GetPair, namely that for every
$(s,a)\in S^{\Diamond}\times{Av}$ we get a converging subsequence
$(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})$ where additionally
$\Delta(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})$ converges to $\Delta(s,a)$ in
total variation. Intuitively, since infinitely many updates occur infinitely
close to $(s,a)$, its limit lower bound $\mathsf{L}_{\infty}(s,a)$ agrees with
the limit of the updates values
$\lim_{k\to\infty}\Delta(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})\langle\mathsf{L}_{\mathsf{t}_{k}}\rangle$.
Since $\mathsf{L}_{\infty}$ satisfies the fixed point equation and is less or
equal to the value function $\mathcal{V}$, we get the result, since
$\mathcal{V}$ is the smallest fixed point.

##  4 Converging Upper Bounds

In this section, we present the second set of assumptions, allowing us to
additionally compute converging _upper_ bounds. With both lower and upper
bounds, we can quantify the progress of the algorithm and, in particular,
terminate the computation once the bounds are sufficiently close. Therefore,
instead of only providing a semi-decision procedure for reachability, this
algorithm is able to determine the maximal reachability probability up to a
given precision. Thus, we obtain the first algorithm able to handle such
general systems with guarantees on its result. We again present our
assumptions together with a discussion of their necessity (Section 4.1), and
then introduce the subsequent algorithm and prove its correctness (Section
4.2). As expected, obtaining this additional information also requires
additional assumptions. On the other hand, quite surprisingly, we can use the
additional information of upper bounds to actually _speed up_ the computation,
as discussed in Section 5.3.

As before, our approach is inspired by algorithms for finite MDP, in this case
by _Bounded Real-Time Dynamic Programming_ (BRTDP) [43, 12]. BRTDP uses the
same update equations as VI, but iterates both lower and upper bounds. A major
contribution of [12] was to solve the long standing open problem of how to
deal with _end components_. These parts of the state space prevent convergence
of the upper bounds by introducing additional fixpoints of Equation 1. We
direct the interested reader to Section A.2 for further details on BRTDP and
insights on the issue of end components. In the uncountable setting, these
issues arise as well alongside several other, related problems, which we
discuss in Section 4.1.2.

###  4.1 Assumptions

The basic assumptions A1–A4 as well as Lipschitz continuity (Assumption C)
remain unchanged. For Maximum Approximation (A2) and Transition Approximation
(A3), we additionally require that we are able to _over_ -approximate the
respective results. The respective assumptions are denoted by A5 and A6,
respectively, and both over-approximations by $\textsc{Approx}_{\geq}$.
Further, we only require a weakened variant of State-Action Sampling, now
called Assumption B.BRTDP instead of Assumption B.VI. Finally, there is the
new Assumption D called Absorption, addressing the aforementioned issue of end
components.

####  4.1.1 B: Weaker Sampling (Asm. B.BRTDP)

We again assume a GetPair oracle, but, perhaps surprisingly, with _weaker_
assumptions. Instead of requiring it to return ?all? actions, we only require
it to yield ?optimal? actions, respective to a given state-action value
function. We first introduce some notation. Intuitively, we want GetPair to
yield actions which are optimal with respect to the _upper bounds_ computed by
the algorithm. However, these upper bounds potentially change after each
update. Thus, assume that $f_{n}\colon S\times{Av}\to[0,1]$ is an arbitrary
sequence of computable, Lipschitz continuous, (point-wise) monotone decreasing
functions, assigning a value to each state-action pair, and set
$\mathcal{F}=(f_{1},f_{2},\dots)$. For each state $s\in S$, set

${Av}_{\mathcal{F}}(s):=\\{a\in{Av}(s)\mid\forall\varepsilon>0.\ \forall
N\in\mathbb{N}.\ \exists n>N.\\\
{\max}_{a^{\prime}\in{Av}(s)}f_{n}(s,a^{\prime})-f_{n}(s,a)<\varepsilon\\},$

i.e. actions that infinitely often achieve values arbitrarily close to the
optimum of $f_{n}$. Let
$S^{\Diamond}_{\mathcal{F}}=\\{last(\varrho)\mid\varrho\in\mathsf{FPaths}_{\mathcal{M},s_{0}}\mathbin{\cap}(S\times{Av}_{\mathcal{F}})^{*}\times
S\\}$ be the set of all states reachable using these optimal
actions.(4)(4)(4)As in Section 3, we simplify the definition of
$S^{\Diamond}_{\mathcal{F}}$ slightly in order to avoid technical details.
Essentially, we require that GetPair samples densely in
$S^{\Diamond}_{\mathcal{F}}\times{Av}_{\mathcal{F}}$.

B.BRTDP: State-Action Sampling

For any $\varepsilon>0$, $\mathcal{F}$ as above, $s\in
S^{\Diamond}_{\mathcal{F}}$, and $a\in{Av}_{\mathcal{F}}(s)$ we have that
GetPair a.s. eventually yields a pair $(s^{\prime},a^{\prime})$ with
$\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))<\varepsilon$ and
$\delta_{TV}(\Delta(s,a),\Delta(s^{\prime},a^{\prime}))<\varepsilon$.

While this new variant may seem much more involved, it is _weaker_ than its
previous variant, since ${Av}_{\mathcal{F}}(s)\subseteq{Av}(s)$ for each $s\in
S$ and thus also $S^{\Diamond}_{\mathcal{F}}\subseteq S^{\Diamond}$. As such,
it also allows for more practical optimizations, which we briefly discuss in
Section 5.3.

####  4.1.2 D: Absorption

We present our most specific assumption. While it is _not needed_ for
correctness, we require it for convergence of the upper bounds to the value
and thus for termination of the algorithm.

D: Absorption

There exists a _known and decidable_ set $R$ (called _sink_) such that
$\mathcal{V}(s)=0$ for all $s\in R$. Moreover, for any $s\in S$ and strategy
$\pi$ we have
$\mathsf{Pr}_{\mathcal{M},s}^{\pi}[\Diamond(T\mathbin{\cup}R)]=1$.

Intuitively, the assumption requires that for all strategies, the system will
eventually reach a target or a goal state; in other words: It is not possible
to avoid both target and sink infinitely long. Variants of this assumption are
used in numerous settings: On MDP, it is similar to the contraction
assumption, e.g. [8, Chp. 4]; in stochastic game theory (a two-player
extension of MDP) it is called _stopping_ , e.g. [16]; and, using terms from
the theory of the stochastic shortest path problem, we require all strategies
to be _proper_ , see e.g. [10].

This assumption already is important in the finite setting: There, Absorption
is equivalent to the absence of _end components_ , which introduce multiple
solutions of Equation 1. Then, a VI algorithm computing upper bounds can be
?stuck? at a greater fixpoint than the value and thus does not converge [12,
22]. _Any_ procedure using value iteration thus either needs to exclude such
cases or detect and treat them. Aside from end components, which are the only
issue in the finite setting, uncountable systems may feature other complex
behaviour, such as Zeno-like approaching the target closer and closer without
reaching it.

Unfortunately, even just detecting these problems already is difficult. For
the mentioned, restricted setting of probabilistic programs, almost sure
termination is $\Pi_{0}^{2}$-complete [36]. Yet, universal termination with
goal set $T\mathbin{\cup}R$ is exactly what we require for Absorption. So,
already on a restricted setting (together with a _given_ guess for $R$), we
cannot decide whether the assumption holds, let alone treat the underlying
problems. Thus, we decide to exclude this issue and delegate treatment to
specialized approaches.

In summary, while this assumption is indeed restrictive, it is the key point
that allows us to obtain convergent upper bounds and thus an anytime
algorithm. As argued above, an assumption of this kind seems to be _necessary_
to obtain such an algorithm in this generality.

######  Remark 4.1.

These problems do not occur when considering _finite horizon_ or _discounted_
properties, which are frequently used in practice. For details on treating
finite horizon objectives, see Section C.1. Discounted reachability with a
factor of $\gamma<1$ is equivalent to normal reachability where at each step
the system moves into a sink state with probability $(1-\gamma)$. Absorption
is trivially satisfied and our methods are directly applicable.

###  4.2 Assumptions Applied: The Convergent Anytime Algorithm

With our assumptions in place, we are ready to present our adaptation of BRTDP
to the uncountable setting. Compared to VI, we now also store upper bounds,
again using Lipschitz-continuity to extrapolate the stored values. In
particular, together with the definitions of Equation 3 we additionally set

$\mathsf{U}(s,a)={\min}_{(s^{\prime},a^{\prime})\in\mathsf{Sampled}}\left(\widehat{\mathsf{U}}(s^{\prime},a^{\prime})+\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))\cdot
C_{\times}\right).$

We also set $\mathsf{U}(s,a)=1$ if either $\mathsf{Sampled}=\emptyset$ or the
above equation would yield $\mathsf{U}(s,a)>1$.

Algorithm 2 The BRTDP algorithm for MDPs with general state- and action-
spaces.

1:ApproxBounds query with precision $\varepsilon$, satisfying A1–A6, B.BRTDP,
C and D.

2:$\varepsilon$-optimal values $(l,u)$.

3:$\mathsf{Sampled}\leftarrow\emptyset,\mathsf{t}\leftarrow 1$
$\triangleright$ Initialize

4:while
$\textsc{Approx}_{\geq}\Big{(}\mathsf{U}(s_{0}),\textsc{Precision}(\mathsf{t})\Big{)}-\textsc{Approx}_{\leq}\Big{(}\mathsf{L}(s_{0}),\textsc{Precision}(\mathsf{t})\Big{)}\geq\varepsilon$
do

5: $s,a\leftarrow\textsc{GetPair}$ $\triangleright$ Sample stat-action pair

6: if $s\in T$ then $\widehat{\mathsf{L}}(s,\cdot)\leftarrow 1$
$\triangleright$ Handle special cases

7: else if $s\in R$ then $\widehat{\mathsf{U}}(s,\cdot)\leftarrow 0$

8: else$\triangleright$ Update upper and lower bounds

9:
$\widehat{\mathsf{U}}(s,a)\leftarrow\textsc{Approx}_{\geq}(\Delta(s,a)\langle\mathsf{U}\rangle,\textsc{Precision}(\mathsf{t}))$

10:
$\widehat{\mathsf{L}}(s,a)\leftarrow\textsc{Approx}_{\leq}(\Delta(s,a)\langle\mathsf{L}\rangle,\textsc{Precision}(\mathsf{t}))$

11: $\mathsf{Sampled}\leftarrow\mathsf{Sampled}\mathbin{\cup}\\{(s,a)\\}$,
$\mathsf{t}\leftarrow\mathsf{t}+1$

12:return $(\mathsf{L}(s_{0}),\mathsf{U}(s_{0}))$

We present BRTDP in Algorithm 2. It is structurally similar to BRTDP in the
finite setting (see Section A.2). The major difference is given by the storage
tables $\widehat{\mathsf{U}}$ and $\widehat{\mathsf{L}}$ used to compute the
current bounds $\mathsf{U}$ and $\mathsf{L}$, again exploiting Lipschitz
continuity. As before, the central idea is to repeatedly update state-action
pairs given GetPair. If GetPair yields a state of the terminal sets $T$ and
$R$, we update the stored values directly. Otherwise, we back-propagate the
value of the selected pair by computing the expected value under this
transition. Moreover, we again require that $\textsc{Precision}(\mathsf{t})$
converges to zero. Note that the algorithm can easily be supplied with
a-priori knowledge by initializing the upper and lower bounds to non-trivial
values. Moreover, in contrast to VI, this algorithm is an _anytime_ algorithm,
i.e. it can at any time provide an approximate solution together with its
precision.

Despite the algorithm being structurally similar to the finite variant of
[12], the proof of correctness unsurprisingly is more intricate due to the
uncountable sets. We again provide both a simplified proof sketch here and the
full technical proof in Section E.2.

######  Theorem 4.2.

Algorithm 2 is correct under Assumptions A1–A6, B.BRTDP, C and D, and
terminates with probability 1.

######  Proof 4.3 (Proof sketch).

We again obtain monotonicity of the bounds, i.e.
$\mathsf{L}_{\mathsf{t}}(s,a)\leq\mathsf{L}_{\mathsf{t}+1}(s,a)\leq\mathcal{V}(s,a)\leq\mathsf{U}_{\mathsf{t}+1}(s,a)\leq\mathsf{U}_{\mathsf{t}}(s,a)$
by induction on $\mathsf{t}$, using completely analogous arguments.

By monotonicity, we also obtain well defined limits $\mathsf{U}_{\infty}$ and
$\mathsf{L}_{\infty}$. Further, we define the difference function
$\mathrm{Diff}_{\mathsf{t}}(s,a)=\mathsf{U}_{\mathsf{t}}(s,a)-\mathsf{L}_{\mathsf{t}}(s,a)$
together with its state based counterpart $\mathrm{Diff}_{\mathsf{t}}(s)$ and
its limit $\mathrm{Diff}_{\infty}(s)$. We show that
$\mathrm{Diff}_{\infty}(s_{0})=0$, proving convergence. To this end, similar
to the previous proof, we prove that $\mathrm{Diff}_{\infty}$ satisfies a
fixed point equation on $S^{\Diamond}_{+}$ (see B.BRTDP), namely
$\mathrm{Diff}_{\infty}(s)=\Delta(s,a(s))\langle\mathrm{Diff}_{\infty}\rangle$
where $a(s)$ is a specially chosen ?optimal? action for each state satisfying
$\mathrm{Diff}_{\infty}(s,a(s))=\mathrm{Diff}_{\infty}(s)$. Now, set
$\mathrm{Diff}_{*}=\max_{s\in S^{\Diamond}_{+}}\mathrm{Diff}_{\infty}(s)$ the
maximal difference on $S^{\Diamond}_{+}$ and let $S^{\Diamond}_{*}$ be the set
of witnesses obtaining $\mathrm{Diff}_{*}$. Then,
$\Delta(s,a(s),S^{\Diamond}_{*})=1$: If a part of the transition’s probability
mass would move to a region with smaller difference, an appropriate update of
a pair close to $(s,a(s))$ would reduce its difference. Hence, the set of
states $S^{\Diamond}_{*}$ is a ?stable? subset of the system when following
the actions $a(s)$. By Absorption, we eventually have to reach either the
target $T$ or the sink $R$ starting from any state in $S^{\Diamond}_{*}$.
Since $\mathrm{Diff}_{\infty}(s)=0$ for all (sampled) states in
$T\mathbin{\cup}R$ and $\mathrm{Diff}_{\infty}$ satisfies the fixed point
equation, we get that $\mathrm{Diff}_{\infty}(s)=0$ for all states
$S^{\Diamond}_{*}$ and consequently $\mathrm{Diff}_{\infty}(s_{0})=0$.

##  5 Discussion

###  5.1 Relation to Algorithms for Finite Systems and Discretization

Our algorithm directly generalizes the classical value iteration as well as
BRTDP for finite MDP by an appropriate choice of GetPair. In value iteration,
it proceeds in round-robin fashion, enumerating all state-action pairs. Note
that the algorithm immediately uses the results of previous updates,
corresponding to the _Gauß-Seidel variant_ of VI; to exactly obtain
synchronous value iteration, we would have to slightly modify the structure
for saving the values. In BRTDP, GetPair simulates paths through the MDP and
we update only those states encountered during the simulation.

Approaches based on discretization through, e.g., grids with increasing
precision, essentially reduce the uncountable state space to a finite one.
This is also encompassed by GetPair, e.g. by selecting the grid points in
round robin or randomized fashion. However, our algorithm has the following
key advantages when compared to classical discretization. Firstly, it avoids
the need to grid the whole state space (typically into cells of regular
sizes). Secondly, in discretization, updating the value of one cell does not
directly affect the value in other cells; in contrast in our algorithm,
knowledge about a state fluently propagates to other areas (by using Equation
3) without being hindered by (arbitrarily chosen) cell boundaries.

###  5.2 Extensions

We outline possible extensions and augmentations of our approach to showcase
its versatility.

Discontinuities Our Lipschitz assumption C actually is slightly stronger than
required. We first give an example of a system exhibiting discontinuities and
then describe how our approach can be modified to deal with it. More details
are in Section C.2.

######  Example 5.1.

Consider a robot navigating a terrain with cliffs, where falling down a cliff
immediately makes it impossible to reach the target. There, states which are
barely on the edge may still reach the goal with significant probability,
while a small step to the side results in falling down the cliff and zero
probability of reaching the goal.

To solve this example, one could model the cliff as a steep but continuous
slope, which would make our approach still possible. Unfortunately, this might
not be very practical, since the Lipschitz constant then is quite large.

However, if we know of discontinuities, e.g. the location of cliffs in the
terrain the robot navigates, both our algorithms can be extended as follows:
Instead of requiring $\mathcal{V}$ to be continuous on the whole domain, we
may assume that we are given a (finite, decidable) partitioning of the state
set $S$ into several sets $S_{i}$. We allow the value function to be
discontinuous along the boundaries of $S_{i}$ (the cliffs), as long as it
remains Lipschitz-continuous inside each $S_{i}$. We only need to slightly
modify the assumption on GetPair by requiring that for any state-action pair
$(s,a)$ with $s\in S_{i}$ we eventually get a nearby, similarly behaving
state-action pair $(s^{\prime},a^{\prime})$ _of the same region_ , i.e.
$s^{\prime}\in S_{i}$. While computing the bounds of a particular state-action
pair, e.g. $\mathsf{U}(s,a)$, we first determine which partition $S_{i}$ the
state $s$ belongs to and then only consider the stored values of states inside
the region $S_{i}$.

Linear Temporal Logic In [12], the authors extend BRTDP to LTL queries [46].
Several difficulties arise in the uncountable setting. For example, in order
to prove _liveness_ conditions, we need to solve the _repeated reachability_
problem, i.e. whether a particular set of states is reached infinitely often.
This is difficult even for restricted classes of uncountable systems, and
impossible in the general case. In particular, [12] relies on analysing end
components, which we already identified as an unresolved problem. We provide
further insight in Section C.3. Nevertheless, there is a straightforward
extension of our approach to the subclass of _reach-avoid_ problems [54] (or
_constrained reachability_ [56]), see Section C.4.

###  5.3 Implementation and Heuristics

For completeness, we implemented a prototype of our BRTDP algorithm to
demonstrate its effectiveness. See Appendix D for details and an evaluation on
both a one- and two-dimensional navigation model. Our implementation is barely
optimized, with no delegation to high-performance libraries. Yet, these non-
trivial models are solved in reasonable time. However, since we aim for
assumptions that are as general as possible, one cannot expect our generic
approach perform on par with highly optimized tools. Our prototype serves as a
proof-of-concept and does not aim to be competitive with specialized
approaches. We highlight again that the goal of our paper is _not_ to be
practically efficient in a particular, restricted setting, but rather to
provide general assumptions and theoretical algorithms applicable to all kinds
of uncountable systems.

Aside from several possible optimizations concerning the concrete
implementation, we suggest two more general directions for heuristics:

Adaptive Lipschitz constants As an example, suppose that a robot is navigating
mostly flat land close to its home, but more hilly terrain further away. The
flat land has a smaller Lipschitz constant than the hilly terrain, and thus
here we can infer tighter bounds. More generally, given a partitioning of the
state space and local Lipschitz constants for every subset, we use this local
knowledge when computing $\widehat{\mathsf{L}}$ and $\widehat{\mathsf{U}}$
instead of using the global Lipschitz constant, which is the maximum of all
local ones. See Section C.2 for details.

GetPair-heuristics In Section 3.1.2, we mentioned two simple implementations
of GetPair. Firstly, we can discretize both state and action space, yielding
each state-action pair in the discretization for a finite number of
iterations, choosing a finer discretization constant, and repeating the
process until convergence. Assuming that we can sample all state-action pairs
in the discretization, this method eventually samples arbitrarily close to
_any_ state-action pair in $S\times{Av}$ and thus trivially satisfies the
sampling assumption. This intuitively corresponds to executing interval
iteration [22] on the (increasingly refined) discretized systems. Note that
this approach completely disregards the reachability probability of certain
states and invests the same computational effort for all of them. In
particular, it invests the same amount of computational effort into regions
which are only reached with probability $10^{-100}$ as in regions around the
initial state $s_{0}$.

Thus, a second approach is to sample a path through the system at random,
following random actions. This approach updates states roughly proportional to
the probability of being reached, which already in the finite setting yields
dramatic speed-ups [38].

However, we can also use further information provided by the algorithm, namely
the upper bounds. As mentioned in [12], following ?promising? actions with a
large upper bound proves to be beneficial, since actions with small upper
bound likely are suboptimal. To extend this idea to the general domain, we
need to apply a bit of care. In particular, it might be difficult to select
exactly from the optimal set of actions, since already
$\operatorname*{arg\,max}_{a\in{Av}(s)}\mathsf{U}(s,a)$ might be very
difficult to compute. Yet, it is sufficient to choose some constant $\xi>0$
and over-approximate the set of $\xi$-optimal actions in a given state,
randomly selecting from this set. This over-approximation can easily be
performed by, for example, randomly sampling the set of available actions
${Av}(s)$ until we encounter an action close to the optimum (which can
approximate due to our assumptions). By generating paths only using these
actions, we combine the previous idea of focussing on ?important? states (in
terms of reachability) with an additional focus on ?promising? states (in
terms of upper bounds). This way, the algorithm _learns_ from its experiences,
using it as a guidance for future explorations.

More generally, we can easily apply more sophisticated learning approaches by
interleaving it with one of the above methods. For example, by following the
learning approach with probability $\nu$ and a ?safe? method with probability
$1-\nu$ we still obtain a safe heuristic, since the assumption only requires
limit behaviour. As such, we can combine our approach with existing, learning
based algorithm by following their suggested heuristic and interleave it with
some sampling runs guided by the above ideas. In other words, this means that
the learning algorithm can focus on finding a reasonable solution quickly,
which is then subsequently _verified_ by our approach, potentially _improving_
the solution in areas where the learner is performing suboptimally. On top,
the (guaranteed) bounds identified by our algorithm can be used as feedback to
the learning algorithm, creating a positive feedback loop, where both
components improve each other’s behaviour and performance.

##  6 Conclusion

In this work, we have presented the first anytime algorithm to tackle the
reachability problem for MDP with uncountable state- and action-spaces, giving
both correctness and termination guarantees under general assumptions. The
experimental evaluation of our prototype implementation shows both promising
results and room for improvements.

On the theoretical side, we conjecture that Assumption D: Absorption can be
weakened if we complement it with an automatic procedure that finds and treats
problematic parts of the state space of a certain kind, similar to the
collapsing approach on finite MDP [22, 12]. Note that as the general problem
is undecidable, some form of Absorption will remain necessary. On the
practical side, we aim for a more sophisticated tool, applying our theoretical
foundation to the full range of MDP, including discrete discontinuities.
Moreover, we want to combine the tool with existing ways of identifying the
Lipschitz constant.

## References

  * [1] Alessandro Abate, Saurabh Amin, Maria Prandini, John Lygeros, and Shankar Sastry.  Computational approaches to reachability analysis of stochastic hybrid systems.  In HSCC, volume 4416 of Lecture Notes in Computer Science, pages 4–17. Springer, 2007.  doi:10.1007/978-3-540-71493-4\\_4.
  * [2] Alessandro Abate, Joost-Pieter Katoen, John Lygeros, and Maria Prandini.  Approximate model checking of stochastic hybrid systems.  Eur. J. Control, 16(6):624–641, 2010.  doi:10.3166/ejc.16.624-641.
  * [3] 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.  doi:10.1016/j.automatica.2008.03.027.
  * [4] Christel Baier and Joost-Pieter Katoen.  Principles of model checking.  MIT Press, 2008.
  * [5] Christel Baier, Joachim Klein, Linda Leuschner, David Parker, and Sascha Wunderlich.  Ensuring the reliability of your model checker: Interval iteration for Markov decision processes.  In CAV (1), volume 10426 of Lecture Notes in Computer Science, pages 160–180. Springer, 2017.
  * [6] Richard Bellman.  Dynamic programming.  Science, 153(3731):34–37, 1966.
  * [7] Dimitri Bertsekas.  Convergence of discretization procedures in dynamic programming.  IEEE Transactions on Automatic Control, 20(3):415–419, 1975.
  * [8] Dimitri P Bertsekas and Steven Shreve.  Stochastic optimal control: the discrete-time case.  1978\.
  * [9] Dimitri P. Bertsekas and John N. Tsitsiklis.  Parallel and distributed computation.  Prentice Hall, 1989.
  * [10] Dimitri P. Bertsekas and John N. Tsitsiklis.  An analysis of stochastic shortest path problems.  Math. Oper. Res., 16(3):580–595, 1991.  doi:10.1287/moor.16.3.580.
  * [11] Patrick Billingsley.  Probability and Measure, volume 939.  John Wiley & Sons, 2012.
  * [12] Tomás Brázdil, Krishnendu Chatterjee, Martin Chmelik, Vojtech Forejt, Jan Kretínský, Marta Z. Kwiatkowska, David Parker, and Mateusz Ujma.  Verification of Markov decision processes using learning algorithms.  In ATVA, volume 8837 of Lecture Notes in Computer Science, pages 98–114. Springer, 2014.  doi:10.1007/978-3-319-11936-6\\_8.
  * [13] John L. Bresina, Richard Dearden, Nicolas Meuleau, Sailesh Ramakrishnan, David E. Smith, and Richard Washington.  Planning under continuous time and resource uncertainty: A challenge for AI.  CoRR, abs/1301.0559, 2013.  URL: http://arxiv.org/abs/1301.0559, arXiv:1301.0559.
  * [14] Debasish Chatterjee, Eugenio Cinquemani, and John Lygeros.  Maximizing the probability of attaining a target prior to extinction.  Nonlinear Analysis: Hybrid Systems, 5(2):367–381, 2011.
  * [15] Krishnendu Chatterjee, Zuzana Kretínská, and Jan Kretínský.  Unifying two views on multiple mean-payoff objectives in Markov decision processes.  Logical Methods in Computer Science, 13(2), 2017.  doi:10.23638/LMCS-13(2:15)2017.
  * [16] Anne Condon.  The complexity of stochastic games.  Inf. Comput., 96(2):203–224, 1992.  doi:10.1016/0890-5401(92)90048-K.
  * [17] Zhengzhu Feng, Richard Dearden, Nicolas Meuleau, and Richard Washington.  Dynamic programming for structured continuous Markov decision problems.  In UAI, pages 154–161. AUAI Press, 2004.
  * [18] Vojtech Forejt, Marta Z. Kwiatkowska, Gethin Norman, and David Parker.  Automated verification techniques for probabilistic systems.  In SFM, volume 6659 of Lecture Notes in Computer Science, pages 53–113. Springer, 2011.  doi:10.1007/978-3-642-21455-4\\_3.
  * [19] Hongfei Fu and Krishnendu Chatterjee.  Termination of nondeterministic probabilistic programs.  In VMCAI, volume 11388 of Lecture Notes in Computer Science, pages 468–490. Springer, 2019.  doi:10.1007/978-3-030-11245-5\\_22.
  * [20] Kush Grover, Jan Kretínský, Tobias Meggendorfer, and Maximilian Weininger.  Anytime guarantees for reachability in uncountable markov decision processes.  CoRR, abs/2008.04824, 2020.  URL: https://arxiv.org/abs/2008.04824.
  * [21] Carlos Guestrin, Milos Hauskrecht, and Branislav Kveton.  Solving factored MDPs with continuous and discrete variables.  In UAI, pages 235–242. AUAI Press, 2004.
  * [22] Serge Haddad and Benjamin Monmege.  Interval iteration algorithm for MDPs and IMDPs.  Theor. Comput. Sci., 735:111–131, 2018.  doi:10.1016/j.tcs.2016.12.003.
  * [23] Sofie Haesaert, Sadegh Soudjani, and Alessandro Abate.  Temporal logic control of general Markov decision processes by approximate policy refinement.  In ADHS, volume 51 of IFAC-PapersOnLine, pages 73–78. Elsevier, 2018.  doi:10.1016/j.ifacol.2018.08.013.
  * [24] Sofie Haesaert, Sadegh Esmaeil Zadeh Soudjani, and Alessandro Abate.  Verification of general Markov decision processes by approximate similarity relations and policy refinement.  SIAM J. Control and Optimization, 55(4):2333–2367, 2017.  doi:10.1137/16M1079397.
  * [25] Arnd Hartmanns and Benjamin Lucien Kaminski.  Optimistic value iteration.  In CAV (2), volume 12225 of Lecture Notes in Computer Science, pages 488–511. Springer, 2020.
  * [26] Mohammadhosein Hasanbeig, Alessandro Abate, and Daniel Kroening.  Certified reinforcement learning with logic guidance.  CoRR, abs/1902.00778, 2019.  URL: http://arxiv.org/abs/1902.00778, arXiv:1902.00778.
  * [27] Mohammadhosein Hasanbeig, Alessandro Abate, and Daniel Kroening.  Logically-constrained neural fitted q-iteration.  In AAMAS, pages 2012–2014. International Foundation for Autonomous Agents and Multiagent Systems, 2019.  URL: http://dl.acm.org/citation.cfm?id=3331994.
  * [28] William B. Haskell, Rahul Jain, Hiteshi Sharma, and Pengqian Yu.  A universal empirical dynamic programming algorithm for continuous state MDPs.  IEEE Trans. Automat. Contr., 65(1):115–129, 2020.  doi:10.1109/TAC.2019.2907414.
  * [29] Thomas A. Henzinger, Peter W. Kopke, Anuj Puri, and Pravin Varaiya.  What’s decidable about hybrid automata?  In STOC, pages 373–382. ACM, 1995.
  * [30] Thomas A. Henzinger, Peter W. Kopke, Anuj Puri, and Pravin Varaiya.  What’s decidable about hybrid automata?  J. Comput. Syst. Sci., 57(1):94–124, 1998.  doi:10.1006/jcss.1998.1581.
  * [31] Onésimo Hernández-Lerma and Jean B Lasserre.  Discrete-time Markov control processes: basic optimality criteria, volume 30.  Springer Science & Business Media, 2012.
  * [32] Ronald A Howard.  Dynamic programming and Markov processes.  1960\.
  * [33] Mingzhang Huang, Hongfei Fu, and Krishnendu Chatterjee.  New approaches for almost-sure termination of probabilistic programs.  In Program. Lang. and Sys., volume 11275 of Lecture Notes in Computer Science, pages 181–201. Springer, 2018.  doi:10.1007/978-3-030-02768-1\\_11.
  * [34] 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.  doi:10.1145/3360555.
  * [35] Manfred Jaeger, Peter Gjøl Jensen, Kim Guldstrand Larsen, Axel Legay, Sean Sedwards, and Jakob Haahr Taankvist.  Teaching stratego to play ball: Optimal synthesis for continuous space MDPs.  In ATVA, volume 11781 of Lecture Notes in Computer Science, pages 81–97. Springer, 2019.  doi:10.1007/978-3-030-31784-3\\_5.
  * [36] Benjamin Lucien Kaminski and Joost-Pieter Katoen.  On the hardness of almost-sure termination.  In MFCS, volume 9234 of Lecture Notes in Computer Science, pages 307–318. Springer, 2015.  doi:10.1007/978-3-662-48057-1\\_24.
  * [37] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  Combinatorica, 4(4):373–396, 1984.  doi:10.1007/BF02579150.
  * [38] Jan Kretínský and Tobias Meggendorfer.  Of cores: A partial-exploration framework for Markov decision processes.  In CONCUR, volume 140 of LIPIcs, pages 5:1–5:17. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019.  doi:10.4230/LIPIcs.CONCUR.2019.5.
  * [39] Ratan Lal and Pavithra Prabhakar.  Bounded verification of reachability of probabilistic hybrid systems.  In QEST, volume 11024 of Lecture Notes in Computer Science, pages 240–256. Springer, 2018.  doi:10.1007/978-3-319-99154-2\\_15.
  * [40] Bernard F Lamond and Abdeslem Boukhtouta.  Water reservoir applications of Markov decision processes.  In Handbook of Markov decision processes, pages 537–558. Springer, 2002.
  * [41] Lihong Li and Michael L. Littman.  Lazy approximation for solving continuous finite-horizon MDPs.  In AAAI, pages 1175–1180. AAAI Press / The MIT Press, 2005\.  URL: http://www.aaai.org/Library/AAAI/2005/aaai05-186.php.
  * [42] Masoud Mahootchi.  Storage system management using reinforcement learning techniques and nonlinear models.  2009\.
  * [43] H. Brendan McMahan, Maxim Likhachev, and Geoffrey J. Gordon.  Bounded real-time dynamic programming: RTDP with monotone upper bounds and performance guarantees.  In ICML, volume 119 of ACM International Conference Proceeding Series, pages 569–576. ACM, 2005.  doi:10.1145/1102351.1102423.
  * [44] Francisco S. Melo, Sean P. Meyn, and M. Isabel Ribeiro.  An analysis of reinforcement learning with function approximation.  In ICML, volume 307 of ACM International Conference Proceeding Series, pages 664–671. ACM, 2008.  doi:10.1145/1390156.1390240.
  * [45] Goran Peskir and Albert Shiryaev.  Optimal stopping and free-boundary problems.  Springer, 2006.
  * [46] Amir Pnueli.  The temporal logic of programs.  In 18th Annual Symposium on Foundations of Computer Science, Providence, Rhode Island, USA, 31 October - 1 November 1977, pages 46–57. IEEE Computer Society, 1977.  doi:10.1109/SFCS.1977.32.
  * [47] Martin L. Puterman.  Markov Decision Processes: Discrete Stochastic Dynamic Programming.  Wiley Series in Probability and Statistics. Wiley, 1994.  doi:10.1002/9780470316887.
  * [48] Tim Quatmann and Joost-Pieter Katoen.  Sound value iteration.  In CAV (1), volume 10981 of Lecture Notes in Computer Science, pages 643–661. Springer, 2018.
  * [49] Scott Sanner, Karina Valdivia Delgado, and Leliane Nunes de Barros.  Symbolic dynamic programming for discrete and continuous state MDPs.  In UAI, pages 643–652. AUAI Press, 2011.
  * [50] Claude Elwood Shannon.  Communication in the presence of noise.  Proceedings of the IRE, 37(1):10–21, 1949.
  * [51] Hiteshi Sharma, Mehdi Jafarnia-Jahromi, and Rahul Jain.  Approximate relative value learning for average-reward continuous state MDPs.  In UAI, page 341. AUAI Press, 2019.  URL: http://auai.org/uai2019/proceedings/papers/341.pdf.
  * [52] Fedor Shmarov and Paolo Zuliani.  Probreach: verified probabilistic delta-reachability for stochastic hybrid systems.  In HSCC, pages 134–139. ACM, 2015.
  * [53] Sadegh Esmaeil Zadeh Soudjani and Alessandro Abate.  Adaptive gridding for abstraction and verification of stochastic hybrid systems.  In QEST, pages 59–68. IEEE Computer Society, 2011.  doi:10.1109/QEST.2011.16.
  * [54] Sean Summers and John Lygeros.  Verification of discrete time stochastic hybrid systems: A stochastic reach-avoid decision problem.  Automatica, 46(12):1951–1961, 2010.  doi:10.1016/j.automatica.2010.08.006.
  * [55] Ilya Tkachev, Alexandru Mereacre, Joost-Pieter Katoen, and Alessandro Abate.  Quantitative automata-based controller synthesis for non-autonomous stochastic hybrid systems.  In HSCC, pages 293–302. ACM, 2013.  doi:10.1145/2461328.2461373.
  * [56] 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.  doi:10.1016/j.ic.2016.11.006.
  * [57] Hado van Hasselt.  Reinforcement learning in continuous state and action spaces.  In Reinforcement Learning, volume 12 of Adaptation, Learning, and Optimization, pages 207–251. Springer, 2012.  doi:10.1007/978-3-642-27645-3\\_7.
  * [58] Moshe Y. Vardi and Pierre Wolper.  Automata-theoretic techniques for modal logics of programs.  J. Comput. Syst. Sci., 32(2):183–221, 1986.  doi:10.1016/0022-0000(86)90026-7.
  * [59] Luis Gustavo Rocha Vianna, Scott Sanner, and Leliane Nunes de Barros.  Continuous real time dynamic programming for discrete and continuous state MDPs.  In 2014 Brazilian Conference on Intelligent Systems, BRACIS 2014, Sao Paulo, Brazil, October 18-22, 2014, pages 134–139. IEEE Computer Society, 2014.  doi:10.1109/BRACIS.2014.34.

##  Appendix A Further Insights

###  A.1 Details on Finite Value Iteration

Value iteration is a technique used to solve, e.g. reachability queries in the
finite MDP setting. It essentially amounts to applying _Bellman iteration_ [6]
(see Equation 2) corresponding to the fixed point equation in Equation 1 [18,
Sec. 4.2]. It is known that on finite MDP this iteration converges to the true
value $\mathcal{V}$ in the limit from below, i.e. for all states $s$ we have
(i) $\lim_{n\to\infty}v_{n}(s)=\mathcal{V}(s)$ and (ii) $v_{n}(s)\leq
v_{n+1}(s)\leq\mathcal{V}(s)$ for all iterations $n$ [47, Thm.
7.2.12](5)(5)(5)Note that reachability is a special case of _expected total
reward_ , obtained by assigning a one-time reward of $1$ to each goal state..
It is not difficult to construct a system where convergence up to a given
precision takes exponential time, but in practice VI often is much faster than
methods based on _linear programming_ (LP) [4, Thm. 10.105], which in theory
has worst-case polynomial runtime and yields precise answers [37]. An
important practical issue of VI is the absence of a _stopping criterion_ ,
i.e. a straightforward way of determining in general whether the current
values $v_{n}(s)$ are close to the true value function $\mathcal{V}(s)$, as
discussed in, e.g. [18, Sec. 4.2].

Traditionally, value iteration computes a value for every state of the system
at each step. However, the iteration can also be executed _asynchronously_.
There, the update order may be chosen by heuristics, as long as fairness
constraints are satisfied, i.e. eventually all states get updated (see e.g.
[9]). Additionally, while VI typically assigns values to each state, one can
also store and update values for each state-action pair separately, i.e.
updating $v_{n+1}(s,a)=\Delta(s,a)\langle v_{n}\rangle$ and derive $v_{n}(s)$
by computing $v_{n}(s):=\max_{a\in{Av}(s)}v_{n}(s,a)$.

###  A.2 Details on BRTDP for Finite MDP

Algorithm 3 The BRTDP algorithm for finite MDP without end components.

1:MDP $\mathcal{M}$, state $s_{0}$, precision $\varepsilon$, target set $T$,
sink set $R$, sampler GetPair.

2:$\varepsilon$-optimal values $(l,u)$.

3:$\mathsf{Sampled}\leftarrow\emptyset$, $s_{1}\leftarrow s_{0}$,
$\mathsf{U}(\cdot,\cdot)\leftarrow 1$, $\mathsf{L}(\cdot,\cdot)\leftarrow 0$
$\triangleright$ Initialize

4:while $\mathsf{U}(s_{0})-\mathsf{L}(s_{0})\geq\varepsilon$ do

5: $s,a\leftarrow\textsc{GetPair}(\mathsf{U})$ $\triangleright$ Get next
state-action pair

6: if $s\in T$ then $\mathsf{L}(s,\cdot)\leftarrow 1$ $\triangleright$ Handle
special cases

7: else if $s\in R$ then $\mathsf{U}(s,\cdot)\leftarrow 0$

8: else$\triangleright$ Update upper and lower bounds

9: $\mathsf{U}(s,a)\leftarrow\Delta(s,a)\langle\mathsf{U}\rangle$

10: $\mathsf{L}(s,a)\leftarrow\Delta(s,a)\langle\mathsf{L}\rangle$

11: $\mathsf{Sampled}\leftarrow\mathsf{Sampled}\mathbin{\cup}\\{(s,a)\\}$

12:return $(\mathsf{L}(s_{0}),\mathsf{U}(s_{0}))$

We briefly summarize the ideas of BRTDP, initially presented in [43] and
further developed in [12]. We present a formal description of BRTDP, adapted
from [12], in Algorithm 3. BRTDP deals with reachability on a _finite_ MDP
$\mathcal{M}=(S,Act,{Av},\Delta)$, i.e. $\lvert S\rvert<\infty$ and $\lvert
Act\rvert<\infty$, with a given target set $T\subseteq S$ and precision
$\varepsilon>0$. The central idea is to apply (asynchronous) value iteration
to compute both lower and upper bounds, iterating until the bounds are
$\varepsilon$-close to each other.

The MDP is assumed to have no _end components_ except in the target set $T$
and a given _sink_ set $R\subseteq S$. Intuitively, end components are parts
of the state space where the system can remain forever under a particular
strategy. For example, suppose there are two states $s_{1},s_{2}$ where
$\Delta(s_{1},a_{1},s_{2})=\Delta(s_{2},a_{2},s_{1})=1$, i.e. the system can
go back and forth between $s_{1}$ and $s_{2}$ indefinitely. Thus
$(\\{s_{1},s_{2}\\},\\{a_{1},a_{2}\\})$ is an end component. Technically, such
end components introduce additional fixed points to the equation of Equation 1
and applying the value iteration Equation 2 to upper bounds would not converge
to the true value function, see [12, Ex. 1] for more details. By excluding end
components (except targets or sinks), we basically get that
$\mathsf{Pr}_{\mathcal{M},s}^{\pi}[\Diamond(T\mathbin{\cup}R)]=1$ for any
state $s$ and strategy $\pi$, i.e. no matter what we do, with probability 1 we
eventually end up in either the sink or the target. The algorithm further
requires that we cannot reach the target once we enter the sink, i.e.
$\mathcal{V}(s)=0$ for all states $s\in R$. Clearly, the upper bounds of all
states in the sink can safely be set to $0$. Since we excluded any additional
fixed point of Equation 1 by our end component assumption, iterations from
below and from above converge to the true value.

The algorithm as presented in [12] repeatedly samples a path until a target or
a sink state is visited and then back-propagates the upper and lower bounds
along this path. In our formulation in Algorithm 3, this is captured by
GetPair as follows: GetPair first samples a path and then returns the states
of the path in reverse order, i.e. starting at the target or sink state and
ending at the initial state. Hence, the updates in 9 and 10 are executed on
all states on the path, and the information of whether a target or sink was
reached is back-propagated. This is why we allow GetPair to be a _stateful_
procedure.

One can also apply this back-propagation globally on all states, which
effectively is done in interval iteration [22]. In contrast to that, sampling
allows the algorithm to focus on ?important? parts of the system, instead of
spending effort on unimportant states. An interesting observation of [12] is
that the upper bounds can guide this sampling efficiently. By choosing actions
with a promising upper bound, we always follow the actions which, given our
current information, could be the best action, a concept sometimes called
_optimism in the face of uncertainty_. We indicate this in Algorithm 3 by
giving GetPair a parameter $\mathsf{U}$.

##  Appendix B Additional Discussion of our Assumptions

###  B.1 Basic Assumptions

We first discuss the standard, usually implicit assumptions.

##### Assumption A1: Metric Space

We require metrics on the state and action spaces in order to define Lipschitz
continuity. More generally, we require a notion of distance to be able to
extrapolate from a particular state to its neighbours. In order to use this
notion in the algorithm, the metric naturally has to be computable. This
assumption is given on practically all reasonable systems, in particular when
considering a well behaved subset of an Euclidean space, e.g. $S\times
Act\subseteq\mathbb{R}^{d}$. Hence, the assumption Metric Space is made
implicitly in many works, with $S$ often assumed to be a subset of
$\mathbb{R}^{d}$. In particular, MDP are often used to model physical
processes, which usually are characterized by real valued variables or, more
generally, variables which allow for a natural notion of distance.

Note that discrete state and action spaces (or state and action spaces with
discrete components) satisfy this assumption, too, using the discrete metric,
i.e. $\operatorname{d}_{S}(x,y)=0$ if $x=y$ and $1$ otherwise. Specifically,
due to our compactness assumption, the discrete parts of the state space
cannot be infinite. So, in essence, this metric lets us investigate each
discrete component separately without extrapolating between them.

##### Assumptions A2, A3, A5, and A6: Maximum Approximation & Transition
Approximation

These assumptions are an immediate consequence of our goal to replicate
Equation 2. Maximum Approximation essentially only requires that we can
compute / approximate ${Av}(s)$ and somehow describe this set. Similarly,
Transition Approximation only imposes some minimal knowledge about
$\Delta(s,a)$. Both assumptions can be realised through, for example, dense
sampling of $f(a)$ and
$g^{\prime}(s^{\prime})\coloneqq\Delta(s,a,s^{\prime})\cdot g(s^{\prime})$.
See Figure 3 for an illustration of the Transition Approximation case. Note
that we only require that $g$ can be approximated up to a certain precision.
However, the overall error introduced by computing these approximations of $g$
can be bounded since we consider probability measures: In particular, suppose
the approximations we obtain are $\hat{g}$, i.e.
$g(x)-\varepsilon\leq\hat{g}(x)\leq g(x)$. Then $\Delta(s,a)\langle
g\rangle-\varepsilon\leq\Delta(s,a)\langle\hat{g}\rangle\leq\Delta(s,a)\langle
g\rangle$. For Maximum Approximation, the ?default? method is even simpler: We
only need to sample actions from ${Av}(s)$ with a distance of at most
$\varepsilon/L_{f}$ (where $L_{f}$ is the Lipschitz constant of $f$). Then, we
can evaluate $f$ at all these positions and know that the true maximum of $f$
is at most $\varepsilon$ larger than the largest sampled value (by Lipschitz
continuity of $f$). The cases of over-approximation are exactly analogous.

$-2$$-1$$0$$1$$2$$0$$5\cdot
10^{-2}$$0.1$$0.15$$s^{\prime}$$\Delta(s,a,s^{\prime})\cdot g(s^{\prime})$
Figure 3:  Illustration to show how we can approximate the successor
expectation $\Delta(s,a)\langle g\rangle$ for a Lipschitz continuous $g$
through sampling. Here, we chose a Gaussian as successor density as and $g(s)$
a zig-zag function. As we know that both functions are Lipschitz continuous,
the product is Lipschitz continuous, too. Hence, the grey area (which can be
computed by evaluating $\Delta(s,a,s^{\prime})\cdot g(s^{\prime})$) gives us a
safe under-approximation of the total integral. By additionally applying the
same idea to $\Delta(s,a,s^{\prime})$ alone, we can safely under-approximate
the total probability mass we have already considered. Together, we obtain a
safe and convergent under-approximation of $\Delta(s,a)\langle g\rangle$. Note
that this only required Lipschitz continuity of the density function of
$\Delta(s,a)$. If we know that we are dealing with, e.g., a Gaussian or a
uniform distribution, we can heavily optimize this process.

##### Assumption A4: Target Computability

Clearly, we need to be able to decide whether a given state is a target state
or not, otherwise the computational problem is not well specified. We
highlight that we do _not_ require an explicit description of the target set
$T$, we only need a procedure to decide $s\in T$. As an example, consider the
uncountable state space of the real numbers in $[0,1]$. Then, a predicate like
$s<\sqrt{(}2)-1$ allows us to decide for every $s\in[0,1]$ whether it is a
target state or not.

###  B.2 Lipschitz Continuity

####  B.2.1 Lipschitz Continuity of Transition Function implies Lipschitz
Continuity of the Value Function

We briefly argue that continuity of the transition function (w.r.t. total
variation) implies continuity of the value function, as similarly shown in,
e.g., [2, Thm. 1] or [53, Thm. 3]. Thus, assume that $\Delta$ is Lipschitz
continuous, i.e. we have that

$\delta_{TV}(\Delta(s,a),\Delta(s^{\prime},a^{\prime}))\leq
L\cdot\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime})).$

Now, recall that

$\mathcal{V}(s,a)\coloneqq\int_{\hat{s}\in
S}\mathcal{V}(\hat{s})\,d\Delta(s,a).$

Since $0\leq\mathcal{V}(s)\leq 1$ for all $s\in S$, we immediately get

$\displaystyle\lVert\mathcal{V}(s,a)-\mathcal{V}(s^{\prime},a^{\prime})\rVert$
$\displaystyle\leq\left\lvert\int_{\hat{s}\in
S}\mathcal{V}(\hat{s})\,d\Delta(s,a)-\int_{\hat{s}\in
S}\mathcal{V}(\hat{s})\,d\Delta(s^{\prime},a^{\prime})\right\rvert$
$\displaystyle\leq\delta_{TV}(\Delta(s,a),\Delta(s^{\prime},a^{\prime}))$
$\displaystyle\leq
L\cdot\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime})).$

To conclude Lipschitz continuity of the value function, one further step is
needed. Recall that
$\mathcal{V}(s)\coloneqq\sup_{a\in{Av}(s)}\mathcal{V}(s,a)$. If the set of
available actions would change abruptly between two nearby states $s$ and
$s^{\prime}$, the continuity of $\mathcal{V}(s,a)$ would not allow us to
conclude anything about the continuity of $\mathcal{V}$.(6)(6)(6)Observe that
this is the underlying reason for our State-Action Sampling assumption, too.
Thus, we furthermore need a ?continuous? behaviour of the action space.
Formally, we require that for each state-action pair $(s,a)$ and state
$s^{\prime}$, there exists an action $a^{\prime}\in{Av}(s^{\prime})$ such that
$\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))\leq
L^{\prime}\cdot\operatorname{d}_{S}(s,s^{\prime})$ for some
$L^{\prime}>0$.(7)(7)(7)This assumption is directly implied by the typical,
much more restrictive assumption of requiring that ${Av}$ is constant. Then,
we can conclude that $\mathcal{V}$ is Lipschitz continuous, too (recall that
$\operatorname{d}_{S}$ and $\operatorname{d}_{\times}$ are compatible).

####  B.2.2 Lipschitz Continuous Value Function

Several works assume Lipschitz continuity of, e.g., the transition function
and derive Lipschitz continuity of the value function, while we only assume
the latter. We demonstrate that our assumption is strictly weaker through a
small example where the transition function is not even continuous, yet the
value function is Lipschitz.

To this end, let $S=Act=[0,1]$, ${Av}(s)=Act$ and $\Delta(s,a)=\\{1\mapsto
1\\}$ if $s<a$, $\\{0\mapsto 1\\}$ if $s>a$, and $\\{s\mapsto 1\\}$ otherwise.
In other words, if $s<a$, we immediately proceed to state $1$, if $s>a$ we go
back to $0$, and for $s=a$ we stay on the spot. As such, a slight change in
either state or action may lead to a large change in the transition dynamics.
More concretely, consider the state $s=0.5$ and the action $a=0.5$. We
slightly change the state by $\varepsilon$ and get
$\Delta(0.5-\varepsilon,0.5)=\\{0\mapsto 1\\}$ or
$\Delta(0.5+\varepsilon,0.5)=\\{1\mapsto 1\\}$, which completely differ from
each other as well as from the original transition
$\Delta(0.5,0.5)=\\{0.5\mapsto 1\\}$. Thus, the transition function is not
continuous.

Yet, when we choose $T=\\{1\\}$, the value function is constant (and hence
Lipschitz continuous), since all states can trivially reach the target by
playing action $1$.

####  B.2.3 Formal Definition of the Frequency Markov Chain

In Example 3.2, we formally consider the following MDP (which actually is a
Markov chain). Let $S=[0,1]\mathbin{\cup}\\{s_{+},s_{-}\\}$ (we add two
distinct state for simplicity, however the example can easily be transformed
to a completely continuous one) and $Act=\\{a\\}$. Now, choose some frequency
$k\in\mathbb{N}$. We define the MDP such that the value function is a periodic
function with frequency $k$. Let $f\colon[0,1]\to[0,1]$ be a Lipschitz
continuous function with $f(0)=f(1)$, e.g., an appropriately scaled sine or
similar. As such, define the transition function such that
$\Delta(s,a)=\\{s_{+}\mapsto f(s\cdot k),s_{-}\mapsto 1-f(s\cdot k)\\}$ for
$0\leq s\leq\frac{1}{k}$ and $\Delta(s,a)=\\{s-\frac{1}{k}\mapsto 1\\}$.
Informally, the states between $0$ and $\frac{1}{k}$ obtain a value according
to $f$ scaled to this interval and all other states simply move $\frac{1}{k}$
to the left. Observe that $\mathcal{V}$ is Lipschitz continuous as long as $f$
is Lipschitz continuous. See Figure 1 for an illustration of the resulting
value function where $f$ is a triangle function and $k=4$. Now, observe that
for any set of finitely many sampled rational points $S^{\prime}\subseteq
S\mathbin{\cap}\mathbb{Q}$, we can choose $k$ and $f$ such that
$\mathcal{V}(s)=1$ for all $s\in S^{\prime}$, but there are also uncountably
many $s^{\prime}$ with $\mathcal{V}(s^{\prime})<\varepsilon$ for every
$\varepsilon>0$. Hence, even though $\mathcal{V}$ is Lipschitz continuous,
without knowing the associated constant we cannot conclude anything about
neighbouring points.

##  Appendix C Extensions and Relaxations

We outline several possible extensions and augmentations of our approach to
showcase the versatility of our assumptions.

###  C.1 Finite horizon (step-bounded) reachability

By using the same idea as in [12], our approach is directly able to handle
finite horizon reachability, also known as _step-bounded_ reachability, i.e.
the probability of reaching a given target set within $n$ steps. We simply
extend all bound functions with a step counter, e.g. $\mathsf{U}(s)$ becomes
$\mathsf{U}(s,i)$, denoting an upper bound on the probability of reaching $T$
within $i$ steps. Similarly, GetPair is supposed to additionally return a step
number $i$ between $0$ and $n$. We then update the lower bound by

$\widehat{\mathsf{L}}(s,a,i)\leftarrow\textsc{Approx}_{\leq}(\Delta(s,a)\langle\mathsf{L}(\cdot,i-1)\rangle,\textsc{Precision}(\mathsf{t}))$

and analogously for the upper bound $\widehat{\mathsf{U}}$. Note that in this
case we do not need the sink set $R$, since after $n$ steps we know that we
will not be able to reach $T$ any more, i.e. $\mathcal{V}(s,n+1)=0$ for all
states $s\in S$.

###  C.2 Discontinuities and Local Lipschitz Constants

0.511.5$0$$0.5$$1$$S$ Figure 4: Example of the function extension on the set
$[0,2]$ when provided with a partitioning $S_{1}=[0,1]$, $S_{2}=[1,1.5]$,
$S_{3}=[1.5,2]$ and different Lipschitz constants $L_{1}=0.5$, $L_{2}=2$,
$L_{3}=1$. As in Figure 2, the dots represent stored values and solid lines
represent the extrapolated functions.

A first, simple extension is given by allowing the value function to be
discontinuous at known locations in a well behaved way. In particular, we
assume that we are given a finite partitioning of the state set $S$ into
several sets $S_{i}$ for $i\in\\{1,\dots,n\\}$, i.e.
$S=\bigcup_{i=1}^{n}S_{i}$ and $S_{i}$ pairwise disjoint. We allow the value
function to be discontinuous along the boundaries of $S_{i}$, as long as it
remains Lipschitz-continuous inside $S_{i}$. Note that we do _not_ require any
special property of $S_{i}$ except decidability and Lipschitz continuity of
the value function in its interior. In particular, $S_{i}$ does not need to be
closed, convex or have a closed form representation. We only need to slightly
modify the assumption on GetPair by requiring that for any state-action pair
$(s,a)$ with $s\in S_{i}$ we eventually get a nearby, similarly behaving
state-action pair $(s^{\prime},a^{\prime})$ where we also have $s^{\prime}\in
S_{i}$.

The only necessary change to the algorithms is the following: While computing
the bounds of a particular state-action pair, e.g. $\mathsf{U}(s,a)$, we first
determine which partition the state belongs to, i.e. find the unique set
$S_{i}$ such that $s\in S_{i}$, and then only consider states sampled inside
this partition. Then, we define

$\mathsf{U}(s,a)={\min}_{(s^{\prime},a^{\prime})\in\mathsf{Sampled},s^{\prime}\in
S_{i}}\left(\widehat{\mathsf{U}}(s^{\prime},a^{\prime})+C_{\times}\cdot\operatorname{d}_{\times}((s,a),(s^{\prime},a^{\prime}))\right)$

and $\mathsf{L}$ analogously, omitting the obvious special cases. Recall that
$\mathsf{Sampled}$ is the set of all state-action pairs sampled so far. It is
easy to see that correctness is preserved. For termination, we can still
establish the respective fixed point equations in the same way.

Similarly, we can allow for local Lipschitz constants for each state, i.e. an
oracle assigns a Lipschitz constant to each state or state-action pair. Then,
whenever we want to compute the bounds of any particular state, we simply use
the individual Lipschitz constants for each state-action pair in
$\mathsf{Sampled}$. To obtain termination, we only require that the Lipschitz
constants are universally bounded.

###  C.3 Linear Temporal Logic

We discuss how _Linear Temporal Logic_ (LTL) [46] can be approached for
uncountable MDP using our approach. We only briefly define LTL and direct the
interested reader to related work, e.g., [4, Chap. 5], [18, Sec. 7.2]. See
[56] for further discussion of LTL on uncountable MDP. Let $AP$ be a finite,
non-empty set of atomic propositions and $a\in AP$ an arbitrary proposition.
Such propositions could, for example, describe ?variable x is larger than 5?
or ?the system is in an unsafe state?. An LTL formula then is given by the
following syntax

$\phi::=a\mid\lnot\phi\mid\phi\land\phi\mid\mathbin{\mathbf{X}}\phi\mid\phi\mathbin{\mathbf{U}}\phi$

with the usual shorthand definitions $\mathsf{false}=a\land\lnot a$,
$\mathsf{true}=\lnot\mathsf{false}$, and $\phi\lor\psi=\lnot(\phi\land\psi)$.
LTL is evaluated over (infinite) sequences of words, i.e. elements of
${(2^{AP})}^{\omega}$. The logical connectives essentially impose restrictions
on the ?current? valuation, i.e. the atomic propositions at the beginning of
the word. The $\mathbin{\mathbf{X}}\phi$ operator requires that $\phi$ holds
in the next step, while $\phi\mathbin{\mathbf{U}}\psi$ demands that the
formula $\phi$ holds at every step until $\psi$ holds. Two common derivations
are $\mathop{\mathbf{F}}\phi:=\mathsf{true}\mathbin{\mathbf{U}}\phi$,
requiring that $\phi$ eventually holds in the future, and
$\mathop{\mathbf{G}}\phi:=\lnot\mathop{\mathbf{F}}\lnot\phi$ requires that
$\phi$ holds globally at every position. As an example,
$a\land(Xb\lor\mathop{\mathbf{F}}\mathop{\mathbf{G}}\lnot a)$ requires that in
the first step we have $a$ and either we have $b$ in the next step, or
eventually we will never see $a$ again, i.e. $a$ only is true finitely often.

In the finite setting, we can equip an MDP with a valuation mapping $\nu\colon
S\to 2^{AP}$, assigning to each state a set of atomic propositions which hold
in this state. This mapping can directly be lifted to paths, i.e. given a path
$\rho$ we derive the respective word
$\nu(\rho)=\nu(\rho(1))\nu(\rho(2))\cdots$. Thus, we can pose _quantitative
LTL queries_ , e.g. ?what is the maximal probability of satisfying the given
formula??. It is known that such queries can be reduced to structural pre-
computations and then solving a reachability query. Based on these ideas, [12]
explains how the BRTDP algorithm can be adapted to accommodate for such an LTL
query.

However, several difficulties arise in the uncountable setting. A central part
of the finite-state algorithm is identifying (winning) end components of the
_product MDP_ (see, e.g., [4, Sec. 10.3, Sec. 10.6.4], for more details on
this _automata-theoretic approach_ [58]), or equivalently solving the
_repeated reachability_ problem. In particular, to prove _safety_ conditions,
for example ?remain inside a region forever?
($\mathop{\mathbf{G}}\texttt{in\\_region}$), one cannot use sampling alone in
general, since even for a single, finite trace it is impossible to give a
positive judgement for such an ?infinite? horizon property. Instead, one needs
to analyse the system’s transition function to infer knowledge about the
infinite horizon behaviour, which is difficult even for restricted classes of
uncountable systems, and impossible in our case, since we treat the transition
function as a black box. See [56, Sec. 4] for further discussion.

Another problem arises already on very simple, ?smooth? systems with equally
simple properties, which we illustrate in the following. Consider an MDP,
where $S=[-2,2]$, ${Av}(s)=[-1,1]$, and
$\Delta(s,a)=\mathrm{unif}([a-1,a+1])$, i.e. uniformly distributed around the
location chosen through $a$. Furthermore, assume that the goal specified by
the LTL formula is to remain in the area $[-1,1]$ forever. Note that this
query is ?stateless?, it is a simple safety requirement. It is easy to see
that by playing action $0$ from every state we satisfy the goal with
probability $1$. Any other strategy which encounters other actions repeatedly
yields an almost sure loss, i.e. the probability of satisfying the goal is
$0$. This particularly shows that obtaining the correct action by sampling has
probability $0$, even though a sampled path following an unsafe strategy may
remain inside the safe area for a very long time. Moreover, the value function
is not continuous, namely $\mathcal{V}(s)=\mathbbm{1}_{[-1,1]}(s)$, since
there is a surely winning strategy for any state in $[-1,1]$. Interestingly,
the state-action value function $\mathcal{V}(s,a)$ is Lipschitz continuous,
namely $\mathcal{V}(s,a)=\frac{1}{2}\min\\{2-a,a+2\\}$ for all $s\in[-1,1]$.
For example, we have that $\mathcal{V}(0,-1)=0.5$, since by playing $-1$ we
only end up in the ?bad? region with $\frac{1}{2}$ probability, otherwise we
can recover by playing optimally.

###  C.4 Reach-avoid problems

Despite that it may seem quite difficult to solve this problem in general even
on simple systems, we actually can apply our approach to so called reach-avoid
problems. These include, for example, a robot navigating towards a recharge
station while avoiding dangerous terrain. More formally, on top of a
reachability query we assume to be given a (measurable and decidable) region
to be avoided $U\subseteq S$.

As in Section 4, we make our usual assumptions, only that we do not require
Lipschitz continuity on the whole state space. Instead, using the ideas of
Section C.2, it is sufficient to assume Lipschitz continuity of the value
function on $S\setminus(T\mathbin{\cup}U)$. Note that we still require the
?sink?-assumptions. For simplicity, assume that
$\mathsf{Pr}_{\mathcal{M},s}^{\pi}[\Diamond(T\mathbin{\cup}U)]=1$ for all
strategies $\pi$ and states $s$. This means that eventually the system either
has to reach the target or will fall into an unrecoverable ?error? state, for
example running out of energy. In this case, our methods are directly
applicable without any major modifications. Note that the central idea is that
(i) we can judge whether a path succeeds or fails based on a finite prefix and
(ii) such a success or failure occurs with probability 1 under any strategy.
We conjecture that our approach is applicable to any system-LTL pair which
satisfies this criterion.

##  Appendix D Evaluation

We implemented a prototype in Java and evaluated it on two models, which
navigate inside one- and two-dimensional state space, respectively. Our
GetPair uses a mixture of global random sampling and path sampling. To reduce
the implementation complexity, our prototype only supports finite action space
and assumes that all states have all actions available. Moreover, currently
the implementation only supports uniform or discrete distributions, however it
easily can be extended to support further types. Under- and over-
approximation is implemented by a (cached) discretized representation of the
lower and upper bound functions. However, adding uncountable actions is not
too different from additional space dimensions. The experiments were carried
out on consumer grade hardware (2.60GHz Intel i7-9750H CPU, 32 GB RAM).

Figure 5:  Experimental results on the one-dimensional model. The bottom part
shows the development of the computed upper and lower bounds over the course
of the algorithm. Furthermore, we show the final result (dashed) as the
average of the computed bounds. The top part shows the development of the
optimal action over time. ?all equal? denotes that all actions are optimal,
?!emergency? that all actions except emergency are optimal, and ?other? refers
to any other combination of optimal actions not explicitly listed.

The first model is an MDP over $S=[-1,1]$ with
$Act=\\{\textsf{idle},\textsf{left},\textsf{right},\textsf{emergency}\\}$. The
target is $T=\\{1\\}$ and the sink is $R=\\{-1\\}$. Essentially, in each state
the system is affected by (i) its own chosen ?thrust? and (ii) the
gravitational forces of the sinks. The latter is computed according to
textbook physics. For the former, the three actions idle, left, and right
yield a small force in the respective direction. The emergency action gives a
significantly higher acceleration towards the right side, however also comes
with the danger of exploding with a probability of $20\%$, realized by moving
to state $-1$ with said probability. The outcome of each action furthermore is
randomized over a continuous interval. Computation terminates after
$236{,}000$ updates, requiring 9 secs and 0.5 GB RAM to achieve the required
precision of $\varepsilon=0.05$ in the initial state $0$.

Figure 5 summarizes the results. When inspecting the optimal actions, four
regions emerge: States close to $-1$ have practically no chance of escaping.
Thus, all actions have a similar value (yellow), close to zero. For the states
around $0$, the most promising action is the emergency thrust (red) to quickly
get outside the gravitational pull of the sink. When slightly above $0$, just
going right (green) without the risk of explosion is safer. In particular, we
identify the trade-off point between emergency and right at approximately
$0.05$, where the bounds exhibit a sharp bend. Lastly, when already close to
$1$, all actions except the emergency are optimal (blue), since we still have
to avoid the emergency thrust to avoid exploding. This includes the left
action, since the thrust of this action is less than the gravitational pull of
the target.

(a) Value Function

(b) Optimal actions

Figure 6:  Experimental results on the two-dimensional model. On the left, we
show the computed value (average of lower and upper bounds). On the right, we
depict the optimal actions, where yellow corresponds to north, green to east,
and purple means that both actions are roughly equivalent.

The second model is conceptually similar to the first one: The system
navigates the state space $S=[0,1]^{2}$ using the actions
$Act=\\{\textsf{north},\textsf{east}\\}$ with the expected outcomes. The
target set is $T=\\{x\mid\lVert x-(1,1)\rVert_{2}<0.05\\}$ with the sink at
$R=\\{x\mid\lVert x-(0.5,0.5)\rVert_{2}<0.05\\}$.(8)(8)(8)We used regions
instead of single points to simplify the implementation of the models. Our
implementation takes 620 seconds, roughly 3 GB, and $605{,}000$ updates to
converge to a precision of $\varepsilon=0.1$.

Figure 6 shows the results for the two-dimensional model. The output is as
expected: Close to the sink, the value is low, as there is little chance of
escaping, and both actions achieve the same values; analogous for the target.
However, in the region around the sink, it is important to choose the right
action in order to navigate around the sink as quickly as possible, see Figure
6b).

##  Appendix E Technical Proofs

##### Definition of Support for General Measures

As mentioned in the preliminaries, $\operatorname{supp}(\mu)$ intuitively
describes the ?smallest? set which $\mu$ assigns a value of $1$. However,
consider, for example, the uniform distribution $\lambda$ over the interval
$[0,1]$. Here, we could always remove single points or even any countable set
(which have a measure of $0$ under $\lambda$) and $\lambda$ still assigns a
value of $1$ to the remaining set. Consequently, the support of $\lambda$
would not be well defined, since from any supposedly ?smallest? set, we could
again remove a single point and obtain an even smaller set. For a well-defined
notion, let $(X,\mathcal{T})$ be a topological space and
$\Sigma_{X}=\mathfrak{B}(\mathcal{T})$ the Borel $\sigma$-algebra. Then, we
define $\operatorname{supp}(\mu)=\\{x\in X\mid x\in
N\in\mathcal{T}\Rightarrow\mu(N)>0\\}$, i.e. the set of points for which any
open neighbourhood $N$ has positive measure. With this definition, we get
$\operatorname{supp}(\lambda)=[0,1]$, as expected.

##### Iteration Counter

In our proofs, we refer to the specific iteration steps of the algorithm and
the values of each variable at these steps. Since every sampled state-action
pair is an atomic step of the algorithm, we use the iteration counter
$\mathsf{t}$ as index of $s$, $a$, $\mathsf{L}$ and $\mathsf{U}$ to refer to
these values. For example, $s_{\mathsf{t}}$ denotes the state sampled in the
$\mathsf{t}$-th step. To ease notation, we assume that any variable which is
not changed in step $\mathsf{t}$ retains its value. For example, any lower
bound $\widehat{\mathsf{L}}_{\mathsf{t}}(s,a)$ of all state-action action
pairs $(s,a)$ which have not been updated in step $\mathsf{t}$ still have the
same value in step $\mathsf{t}+1$. We keep this convention throughout all the
proofs in Appendix E.

###  E.1 Detailed Proof of Theorem 3.3

We first provide a proof sketch, illustrating the main steps of the proof. The
sketch particularly omits several subtle technicalities. For example, the
precision of values obtained by $\textsc{Approx}_{\leq}$ can be non-monotonic
and approximation errors might accumulate. These issues are treated in the
main proof.

######  Proof E.1 (Proof sketch).

First, we show that
$\mathsf{L}_{\mathsf{t}}(s)\leq\mathsf{L}_{\mathsf{t}+1}(s)\leq\mathcal{V}(s)$
by simple induction on the step. Initially, we have $\mathsf{L}_{1}(s)=0$,
obviously satisfying the condition. The updates in Lines 6 and 7 both keep
correctness, i.e. $\mathsf{L}_{\mathsf{t}+1}(s)\leq\mathcal{V}(s)$, proving
the claim.

Since $\mathsf{L}_{\mathsf{t}}$ is monotone as argued above, its limit for
$\mathsf{t}\to\infty$ is well defined, denoted by $\mathsf{L}_{\infty}$. By
State-Action Sampling, the set of accumulation points of $s_{\mathsf{t}}$
contains all reachable states $S^{\Diamond}$. We then prove that
$\mathsf{L}_{\infty}$ satisfies the fixed point equation Equation 1. For this,
we use the second part of the assumption on GetPair, namely that for every
$(s,a)\in S^{\Diamond}\times{Av}$ we get a converging subsequence
$(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})$ where additionally
$\Delta(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})$ converges to $\Delta(s,a)$ in
total variation. Intuitively, since infinitely many updates occur infinitely
close to $(s,a)$, its limit lower bound $\mathsf{L}_{\infty}(s,a)$ agrees with
the limit of the updates values
$\lim_{k\to\infty}\Delta(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})\langle\mathsf{L}_{\mathsf{t}_{k}}\rangle$.
Since $\mathsf{L}_{\infty}$ satisfies the fixed point equation and is less or
equal to the value function $\mathcal{V}$, we get the result, since
$\mathcal{V}$ is the smallest fixed point.

Before we begin with the technical proof of the correctness theorem, we point
out a subtle issue in the update computation of the algorithm. In the proof,
we want to show that the lower bounds are monotonically increasing, i.e.
$\mathsf{L}_{\mathsf{t}}(s,a)\leq\mathsf{L}_{\mathsf{t}+1}(s,a)$ for any
state-action pair $(s,a)\in S\times{Av}$. To obtain this monotonicity in the
presence of under-approximations, we need to slightly modify the algorithm in
a special case. In particular, suppose that we sampled some state action pair
$(s,a)$ in step $\mathsf{t}$ and its lower bound
$\widehat{\mathsf{L}}_{\mathsf{t}+1}$ was approximated very precisely. If we
now again sample $(s,a)$ exactly in $\mathsf{t}+1$, the approximation
computation may yield a coarser result, consequently we would have
$\mathsf{L}_{\mathsf{t}+2}(s,a)<\mathsf{L}_{\mathsf{t}+1}$. To remedy this
situation, we modify the algorithm to update $\widehat{\mathsf{L}}$ with the
maximum of the current and the computed value. Note that if we sample a nearby
state-action pair $(s^{\prime},a^{\prime})$ in the second step, we do not need
this special treatment, since then we would simply have that
$\widehat{\mathsf{L}}_{\mathsf{t}+2}(s^{\prime},a^{\prime})<\mathsf{L}_{\mathsf{t}+2}(s^{\prime},a^{\prime})$,
since $\mathsf{L}_{\mathsf{t}+2}(s^{\prime},a^{\prime})$ is then computed
based on the stored value of the nearby pair $(s,a)$.

######  Proof E.2.

We show that (i) $\mathsf{L}_{\mathsf{t}}(s)\leq\mathcal{V}(s)$ for all $s\in
S$ and steps $\mathsf{t}$ and (ii)
$\lim_{\mathsf{t}\to\infty}\mathsf{L}_{\mathsf{t}}(s_{0})=\mathcal{V}(s_{0})$.

First, we prove that
$\mathsf{L}_{\mathsf{t}}(s)\leq\mathsf{L}_{\mathsf{t}+1}(s)\leq\mathcal{V}(s)$
for all $s\in S$ and steps $\mathsf{t}$ by induction on $\mathsf{t}$.
Initially, we only need to show that
$\mathsf{L}_{1}(s,a)\leq\mathcal{V}(s,a)$. This clearly is the case since we
have $\mathsf{L}_{1}(s,a)=0$, since $\mathsf{Sampled}$ is empty and by
convention we then set the lower bound to 0. Assume we have
$\mathsf{L}_{\mathsf{t}}(s)\leq\mathcal{V}(s)$ for some step $\mathsf{t}$. The
update in Line 6 is obviously correct, since $\mathcal{V}(s)=1$ for all $s\in
T$. Note that in this case we trivially get monotonicity of the updates. For
the back-propagation in Line 7, observe that
$\Delta(s,a)\langle\mathsf{L}_{\mathsf{t}}\rangle\leq\Delta(s,a)\langle\mathcal{V}\rangle$
by induction. Moreover,
$\mathcal{V}(s,a)=\Delta(s,a)\langle\mathcal{V}\rangle$ by definition. Hence
$\widehat{\mathsf{L}}_{\mathsf{t}+1}(s_{\mathsf{t}},a_{\mathsf{t}})\leq\mathcal{V}(s,a)$.
Thus, correctness is preserved. Monotonicity is directly obtained due to the
above discussion.

Now, we prove that if the algorithm does not terminate we have that
$\lim_{\mathsf{t}\to\infty}\mathsf{L}_{\mathsf{t}}(s_{0})=\mathcal{V}(s_{0})$.
Note that when the algorithm does terminate, there is nothing left to prove.
In the following, we only argue using limit behaviour. Hence, because (i) our
definition of $\textsc{Precision}(\mathsf{t})$ ensures that it converges to 0
in the limit and (ii) by Assumptions A2 and A3 we can approximate arbitrarily
precisely, we have that any arising approximation computation is arbitrarily
precise (this in particular means that accumulation of approximation bounds is
not an issue). We first set up some auxiliary notation. Set
$\mathsf{L}_{\infty}(s,a)=\lim_{\mathsf{t}\to\infty}\mathsf{L}_{\mathsf{t}}(s,a)$,
$\mathsf{L}_{\infty}(s)=\max_{a\in{Av}(s)}\mathsf{L}_{\infty}(s,a)$ (note that
both $\mathsf{L}_{\infty}$ are continuous). These limits are well-defined due
to the above result – the functions are bounded and monotone. Moreover, since
the set of state-action pairs is compact, the convergence of
$\mathsf{L}_{\mathsf{t}}$ is uniform. Let further
$S_{\infty}=\\{s\mid\forall\varepsilon>0.\ \forall\mathsf{t}.\
\exists\mathsf{t}^{\prime}>\mathsf{t}.\
\operatorname{d}_{S}(s_{\mathsf{t}^{\prime}},s)<\varepsilon\\}$ the set of all
accumulation points of $s_{\mathsf{t}}$, i.e. all states to which the
algorithm gets arbitrarily close infinitely often. By compactness of $S$, this
set is not empty. Note that despite the set of all sampled states being
countable, $S_{\infty}$ may be uncountable. Next, for each $s\in S_{\infty}$,
set $Act_{\infty}(s)=\\{a\in{Av}(s)\mid\forall\varepsilon>0.\
\forall\mathsf{t}.\ \exists\mathsf{t}^{\prime}>\mathsf{t}.\
\operatorname{d}_{Act}(a,a_{\mathsf{t}^{\prime}})<\varepsilon\\}$. By our
assumption State-Action Sampling, we have that (i) $S_{\infty}$ contains all
reachable states, i.e. $S^{\Diamond}\subseteq S_{\infty}$, and (ii)
$Act_{\infty}(s)={Av}(s)$ for all $s\in S_{\infty}$. We now prove that for all
$s\in S_{\infty}$ we have that either $s\in T$ and $\mathsf{L}_{\infty}(s)=1$
or $\mathsf{L}_{\infty}(s,a)=\Delta(s,a)\langle\mathsf{L}_{\infty}\rangle$.
Note that the first case is trivial by the update rule of the algorithm. Let
$s\in S^{\Diamond}\setminus T$ arbitrary and $a\in{Av}(s)$. From our
assumption, we obtain a sequence of state-action pairs
$(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})$ which converges to $(s,a)$ and
$\Delta(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})$ converges in total variation
to $\Delta(s,a)$. By definition of the algorithm, we have that
$\mathsf{L}_{\mathsf{t}_{k}+1}(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})=\Delta(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})\langle\mathsf{L}_{\mathsf{t}_{k}}\rangle$.
By uniform convergence and continuity of $\mathsf{L}_{\mathsf{t}}$, we get
$\lim_{k\to\infty}\mathsf{L}_{\mathsf{t}_{k}+1}(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})=\lim_{k\to\infty}\mathsf{L}_{\infty}(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})=\mathsf{L}_{\infty}(s,a)$.
By total variation convergence, we obtain that
$\lim_{k\to\infty}\Delta(s_{\mathsf{t}_{k}},a_{\mathsf{t}_{k}})\langle\mathsf{L}_{\infty}\rangle=\Delta(s,a)\langle\mathsf{L}_{\infty}\rangle$.
The desired claim follows.

To conclude, observe that $\mathcal{V}$ is the least fixed point of the
equation system in Equation 1. Since we have shown that
$\mathsf{L}_{\infty}(s)\leq\mathcal{V}(s)$ for all $s\in S^{\Diamond}$ and
that $S^{\Diamond}\subseteq S^{\Diamond}_{*}$, $\mathsf{L}_{\infty}$ equals
the least fixed point on these states.

###  E.2 Proof of Theorem 4.2

In the algorithm, the $f_{n}$ of Assumption B.BRTDP corresponds to the
computed upper bounds at step $n$ and we write $S^{\Diamond}_{+}$ in the
following to denote this set of states generated by the algorithm to simplify
notation.

######  Proof E.3 (Proof sketch).

We again obtain monotonicity of the bounds, i.e.
$\mathsf{L}_{\mathsf{t}}(s,a)\leq\mathsf{L}_{\mathsf{t}+1}(s,a)\leq\mathcal{V}(s,a)\leq\mathsf{U}_{\mathsf{t}+1}(s,a)\leq\mathsf{U}_{\mathsf{t}}(s,a)$
by induction on $\mathsf{t}$, using completely analogous arguments.

By monotonicity, we also obtain well defined limits $\mathsf{U}_{\infty}$ and
$\mathsf{L}_{\infty}$. Further, we define the difference function
$\mathrm{Diff}_{\mathsf{t}}(s,a)=\mathsf{U}_{\mathsf{t}}(s,a)-\mathsf{L}_{\mathsf{t}}(s,a)$
together with its state based counterpart $\mathrm{Diff}_{\mathsf{t}}(s)$ and
its limit $\mathrm{Diff}_{\infty}(s)$. We show that
$\mathrm{Diff}_{\infty}(s_{0})=0$, proving convergence. To this end, similar
to the previous proof, we prove that $\mathrm{Diff}_{\infty}$ satisfies a
fixed point equation on $S^{\Diamond}_{+}$ (see B.BRTDP), namely
$\mathrm{Diff}_{\infty}(s)=\Delta(s,a(s))\langle\mathrm{Diff}_{\infty}\rangle$
where $a(s)$ is a specially chosen ?optimal? action for each state satisfying
$\mathrm{Diff}_{\infty}(s,a(s))=\mathrm{Diff}_{\infty}(s)$. Now, set
$\mathrm{Diff}_{*}=\max_{s\in S^{\Diamond}_{+}}\mathrm{Diff}_{\infty}(s)$ the
maximal difference on $S^{\Diamond}_{+}$ and let $S^{\Diamond}_{*}$ be the set
of witnesses obtaining $\mathrm{Diff}_{*}$. Then,
$\Delta(s,a(s),S^{\Diamond}_{*})=1$: If a part of the transition’s probability
mass would move to a region with smaller difference, an appropriate update of
a pair close to $(s,a(s))$ would reduce its difference. Hence, the set of
states $S^{\Diamond}_{*}$ is a ?stable? subset of the system when following
the actions $a(s)$. By Absorption, we eventually have to reach either the
target $T$ or the sink $R$ starting from any state in $S^{\Diamond}_{*}$.
Since $\mathrm{Diff}_{\infty}(s)=0$ for all (sampled) states in
$T\mathbin{\cup}R$ and $\mathrm{Diff}_{\infty}$ satisfies the fixed point
equation, we get that $\mathrm{Diff}_{\infty}(s)=0$ for all states
$S^{\Diamond}_{*}$ and consequently $\mathrm{Diff}_{\infty}(s_{0})=0$.

Before the full technical proof, we first prove a small auxiliary lemma,
showing that the upper and lower bounds are monotone and sound, i.e.
$\mathsf{L}$ is increasing, $\mathsf{U}$ is decreasing, and the value always
lies between them. Note that we again use the same adaptations as discussed in
Section E.1 in order to ensure monotonicity.

######  Lemma E.4.

For any step $\mathsf{t}$, state $s$, and action $a\in{Av}(s)$, we have that
$\mathsf{L}_{\mathsf{t}}(s,a)\leq\mathsf{L}_{\mathsf{t}+1}(s,a)$,
$\mathsf{U}_{\mathsf{t}+1}(s,a)\leq\mathsf{U}_{\mathsf{t}}(s,a)$, and
$\mathsf{L}_{\mathsf{t}}(s,a)\leq\mathcal{V}(s,a)\leq\mathsf{U}_{\mathsf{t}}(s,a)$.

######  Proof E.5.

We prove by induction on the step $\mathsf{t}$. Initially, we only need to
show that $\mathsf{L}_{1}(s,a)\leq\mathcal{V}(s,a)\leq\mathsf{U}_{1}(s,a)$.
This clearly is the case since we have $\mathsf{L}_{1}(s,a)=0$ and
$\mathsf{U}_{1}(s,a)=1$.

Now, assume we have
$\mathsf{L}_{\mathsf{t}}(s)\leq\mathcal{V}(s)\leq\mathsf{U}_{\mathsf{t}}(s)$
for some step $\mathsf{t}$. The update in Line 6 is obviously correct, since
$\mathcal{V}(s)=1$ for all $s\in T$. The correctness of Line 7 follows
directly from our assumption on $R$ ($\mathcal{V}(s)=0$ for all $s\in R$).
Note that in these two cases we trivially get monotonicity of the updates. For
the back-propagation in Lines 9 and 10, observe that
$\Delta(s,a)\langle\mathsf{L}_{\mathsf{t}}\rangle\leq\Delta(s,a)\langle\mathcal{V}\rangle\leq\Delta(s,a)\langle\mathsf{U}_{\mathsf{t}}\rangle$
by induction. Moreover,
$\mathcal{V}(s,a)=\Delta(s,a)\langle\mathcal{V}\rangle$ by definition. Hence
$\widehat{\mathsf{L}}_{\mathsf{t}+1}(s_{\mathsf{t}},a_{\mathsf{t}})\leq\mathcal{V}(s,a)\leq\widehat{\mathsf{U}}_{\mathsf{t}+1}(s_{\mathsf{t}},a_{\mathsf{t}})$.
Thus, correctness is preserved.

With this lemma, we can now prove correctness and termination of our
algorithm. Correctness follows directly from the above lemma. In order to
prove termination, we essentially construct a contradiction based on the
Absorption assumption.

######  Proof E.6 (Proof of Theorem 4.2).

Correctness: Follows directly from Lemma E.4. In particular, when the
algorithm terminates, we know that
$\mathsf{L}_{\mathsf{t}}(s_{0})\leq\mathcal{V}(s_{0})\leq\mathsf{U}_{\mathsf{t}}(s_{0})$.

Termination: We prove by contradiction. Thus, assume that the algorithm does
not converge, i.e. we have
$\mathsf{U}_{\mathsf{t}}(s_{0})-\mathsf{L}_{\mathsf{t}}(s_{0})>\varepsilon$
for all steps $\mathsf{t}$. As before, we only argue using limit behaviour,
and again assume that any arising approximation computation is arbitrarily
precise. First, we need to set up some auxiliary notation. Set
$\mathsf{U}_{\infty}(s,a)=\lim_{\mathsf{t}\to\infty}\mathsf{U}_{\mathsf{t}}(s,a)$,
$\mathsf{U}_{\infty}(s)=\max_{a\in{Av}(s)}\mathsf{U}_{\infty}(s,a)$ (note that
both $\mathsf{U}_{\infty}$ are continuous), and analogously define
$\mathsf{L}_{\infty}$. All limits are well-defined due to Lemma E.4 – the
functions are bounded and monotone. Moreover, the convergence is uniform due
to pointwise convergence on a compact domain. Now, define the difference
function
$\mathrm{Diff}_{\mathsf{t}}(s,a)=\mathsf{U}_{\mathsf{t}}(s,a)-\mathsf{L}_{\mathsf{t}}(s,a)$.
For any step $\mathsf{t}$ and state $s$, let
$\mathsf{MaxA}_{\mathsf{t}}(s):=\operatorname*{arg\,max}_{a\in{Av}(s)}\mathsf{U}_{\mathsf{t}}(s,a)$
denote the set of $\mathsf{U}$-optimal actions in state $s$ and set
$\mathrm{Diff}_{\mathsf{t}}(s)=\max_{a\in\mathsf{MaxA}_{\mathsf{t}}(s)}\mathrm{Diff}_{\mathsf{t}}(s,a)$
the maximal difference among them. It is important to note that
$\mathrm{Diff}_{\mathsf{t}}(s)$ does not necessarily equal
$\max_{a\in{Av}(s)}\mathrm{Diff}_{\mathsf{t}}(s,a)$, but it is easy to show
that
$\mathrm{Diff}_{\mathsf{t}}(s)\geq\mathsf{U}_{\mathsf{t}}(s)-\mathsf{L}_{\mathsf{t}}(s)$.
Clearly,
$\mathrm{Diff}_{\infty}(s)=\limsup_{\mathsf{t}\to\infty}\mathrm{Diff}_{\mathsf{t}}(s)$
is well defined, too, since $\mathrm{Diff}_{\mathsf{t}}$ is bounded. Now, our
overall proof strategy is to prove that $\mathrm{Diff}_{\infty}(s_{0})=0$
almost surely, since this implies that eventually
$\mathsf{U}_{\mathsf{t}}(s_{0})-\mathsf{L}_{\mathsf{t}}(s_{0})\leq\mathrm{Diff}_{\mathsf{t}}(s)<\varepsilon$,
contradicting our initial assumption.

Let now $\mathsf{MaxA}_{\infty}(s)=\\{a\in{Av}(s)\mid\forall\varepsilon>0.\
\forall\mathsf{t}.\ \exists\mathsf{t}^{\prime}>\mathsf{t}.\
\mathsf{U}_{\mathsf{t}^{\prime}}(s)-\mathsf{U}_{\mathsf{t}^{\prime}}(s,a)<\varepsilon\\}$
all actions which infinitely often achieve a value arbitrarily close to the
optimum. We now show that for any fixed state $s$ there exists an action
$a_{\infty}^{\max}(s)\in\mathsf{MaxA}_{\infty}(s)\subseteq{Av}(s)$ such that
$\mathrm{Diff}_{\infty}(s)=\mathrm{Diff}_{\infty}(s,a)$. For every
$\mathsf{t}$, choose an arbitrary action
$a_{\mathsf{t}}^{\max}\in\operatorname*{arg\,max}_{a\in\mathsf{MaxA}_{\mathsf{t}}(s)}\mathrm{Diff}_{\mathsf{t}}(s,a)$.
We have that
$\mathsf{U}_{\mathsf{t}}(s,a_{\mathsf{t}}^{\max})=\mathsf{U}_{\mathsf{t}}(s)$
by definition. Consequently, the limit
$\lim_{\mathsf{t}\to\infty}\mathsf{U}_{\mathsf{t}}(s,a_{\mathsf{t}}^{\max})$
is well defined and equals $\mathsf{U}_{\infty}(s)$. Moreover,
$\liminf_{\mathsf{t}\to\infty}\mathsf{L}_{\mathsf{t}}(s,a_{\mathsf{t}}^{\max})$
is well defined and we can choose a subsequence of $a_{\mathsf{t}}^{\max}$
obtaining this limit. Also, since ${Av}(s)$ is compact, there exists an
accumulation point $a_{\infty}^{\max}(s)\in{Av}(s)$ of this subsequence. Note
that $a_{\infty}^{\max}(s)\in\mathsf{MaxA}_{\infty}(s)$. Together, we get that

$\displaystyle\mathrm{Diff}_{\infty}(s)$
$\displaystyle=\limsup_{\mathsf{t}\to\infty}\left(\mathsf{U}_{\mathsf{t}}(s,a_{\mathsf{t}}^{\max})-\mathsf{L}_{\mathsf{t}}(s,a_{\mathsf{t}}^{\max})\right)$
$\displaystyle=\lim_{\mathsf{t}\to\infty}\mathsf{U}_{\mathsf{t}}(s,a_{\mathsf{t}}^{\max})-\liminf_{\mathsf{t}\to\infty}\mathsf{L}_{\mathsf{t}}(s,a_{\mathsf{t}}^{\max})$
$\displaystyle=\mathsf{U}_{\infty}(s,a_{\infty}^{\max}(s))-\mathsf{L}_{\infty}(s,a_{\infty}^{\max}(s))$
$\displaystyle=\mathrm{Diff}_{\infty}(s,a_{\infty}^{\max}(s)).$

Next, we will show that for a particular subset of states, we have that
$\mathrm{Diff}$ satisfies the fixed point equation, i.e.
$\mathrm{Diff}_{\infty}(s)=\Delta(s,a_{\infty}^{\max}(s))\langle\mathrm{Diff}_{\infty}\rangle$,
using our assumptions on GetPair. Let thus
$S_{\infty}=\\{s\mid\forall\varepsilon>0.\ \forall\mathsf{t}.\
\exists\mathsf{t}^{\prime}>\mathsf{t}.\
\operatorname{d}_{S}(s,s_{\mathsf{t}^{\prime}})<\varepsilon\\}$ the set of all
accumulation points of $s_{\mathsf{t}}$, i.e. all states to which the
algorithm gets arbitrarily close infinitely often. By compactness of $S$, this