Quantitative Approximations for Reach-Avoid over General Markov Processes

Sampling-based Approximations with Quantitative
Performance for the Probabilistic Reach-Avoid Problem over
General Markov Processes

arXiv:1409.0553v2 [cs.SY] 9 Sep 2015

Sofie Haesaert

s.haesaert@tue.nl

Control Systems Group, Department of Electrical Engineering
Eindhoven University of Technology
Eindhoven, The Netherlands

Robert Babuska

r.babuska@tudelft.nl

Delft Center of Systems and Control
Delft University of Technology
Delft, The Netherlands

Alessandro Abate

alessandro.abate@cs.ox.ac.uk

Department of Computer Science
University of Oxford
Oxford, United Kingdom

Editor:

Abstract
This article deals with stochastic processes endowed with the Markov (memoryless) property and evolving over general (uncountable) state spaces. The models further depend on a
non-deterministic quantity in the form of a control input, which can be selected to affect the
probabilistic dynamics. We address the computation of maximal reach-avoid specifications,
together with the synthesis of the corresponding optimal controllers. The reach-avoid specification deals with assessing the likelihood that any finite-horizon trajectory of the model
enters a given goal set, while avoiding a given set of undesired states. This article newly
provides an approximate computational scheme for the reach-avoid specification based on
the Fitted Value Iteration algorithm, which hinges on random sample extractions, and gives
a-priori computable formal probabilistic bounds on the error made by the approximation
algorithm: as such, the output of the numerical scheme is quantitatively assessed and thus
meaningful for safety-critical applications. Furthermore, we provide tighter probabilistic error bounds that are sample-based. The overall computational scheme is put in relationship
with alternative approximation algorithms in the literature, and finally its performance is
practically assessed over a benchmark case study.
Keywords: General state-space processes, reach-avoid problem, dynamic programming,
fitted value iteration, computational approximation with error bounds.

1. Introduction
This contribution concerns a problem grounded in concepts from a few different areas: we
deal with probabilistic processes evolving over continuous (and in particular uncountable)
state spaces – this leads to the use of measure-theoretical material from Stochastic Processes
and Probability Theory (Meyn and Tweedie, 1993); we work with models endowed with a
1

Haesaert, Babuska and Abate

control input and investigate control synthesis, which relate to a broad literature in Control Theory (Bertsekas and Shreve, 1996); furthermore, we are interested in quantifying the
probability associated to a dynamical property, known as reach-avoid, which corresponds
to a widely used model specification in the field of Formal Verification (Baier and Katoen,
2008); and finally we employ a sampling-based algorithm to approximately compute the
likelihood associated to the above specification. The algorithm, known as Fitted Value Iteration (FVI) (Munos and Szepesvari, 2008), is a regression scheme developed in Machine
Learning.
We focus on stochastic processes endowed with the Markov property (where the future
is independent of the past, conditional on the present) and, aiming for generality, we deal
with processes evolving over a continuous state space. We are further interested in a class
of such models known as stochastic hybrid systems (SHS) (Abate et al., 2008), which are
endowed with a “hybrid” (that is, both continuous and discrete) state space, which are relevant for a number of applications in Engineering and the Life Sciences (Blom and Lygeros,
2006; Cassandras and J. Lygeros, 2006). This work investigates the problem of controller
synthesis over these models, namely the selection of sequences of control inputs (which in
particular can be functions of the states of the model) over a finite time horizon, in order
to optimise a given figure of merit.
As for the figure of merit of interest in this work, we choose to go beyond the classical
properties investigated in Systems and Control theory, which by and large deal with known
and standard problems of stability, regulation, and tracking. Instead, we focus on the reachavoid specification, a property that is well known and central within the Formal Verification
field (Baier and Katoen, 2008). Notice that classical results in Formal Verification deal with
simple models – usually finite-state transition systems or Markov chains – which allow for
the development of computational results, and which mostly deal with verification tasks that
do not involve policy synthesis. In this work instead we consider reach-avoid specifications
over models with continuous stochastic transitions and endowed with control inputs.
The reach-avoid problem deals with computing the likelihood that, within a given finite
time horizon, a trajectory of the model enters a goal set, while avoiding a given set of
undesired states (both sets are arbitrary measurable subsets of the state space). Equivalently, the property can be expressed as the probability that the process enters a goal set
while dwelling within a set of allowed states. The reach-avoid property is a generalization
of widely studied properties, such as reachability and invariance, and represents a known
specification (denoted as “bounded until”) that lies at the core of a number of modal logics
used in the field of formal verification, such as Linear Temporal Logic and Computational
Time Logic (Baier and Katoen, 2008). From a controller synthesis perspective, the goal
becomes that of either maximizing or minimizing the above specification over the given
time horizon.
In the context of probabilistic models evolving over continuous domains and in discrete
time (which is the framework considered in this work), the probabilistic reachability and
reach-avoid specifications have been investigated in (Abate et al., 2008; Summers and Lygeros,
2010). These results have recently led to the study of other properties, either richer
(Abate et al., 2011) or defined over unbounded time horizons (Tkachev and Abate, 2014).
These results have focused on the theoretical characterization of the specifications/properties
2

Quantitative Approximations for Reach-Avoid over General Markov Processes

of interest: of course it is also of much interest to provide algorithms that can numerically
compute these figures. Computational approaches to probabilistic reachability have been
studied in (Abate et al., 2010; Esmaeil Zadeh Soudjani and Abate, 2013): the strength of
these results is that the proposed numerical schemes have explicitly quantified error bounds.
This is unlike other, known approximation schemes in the literature (Koutsoukos and Riley,
2006; Kushner and Dupuis, 2001; Prandini and Hu, 2006), which provide results with properties that are only known asymptotically.
This article provides a new approximate computational scheme for the reach-avoid specification based on the Fitted Value Iteration algorithm, which hinges on random sample
extractions. This work originally derives formal probabilistic bounds on the error made by
the approximation algorithm. In order to do so, the FVI scheme is tailored to the characterization of the reach-avoid problem, which leads to Dynamic Programming (DP) recursions
based on a sum-multiplicative form that is non-standard since it departs from the classical
additive (possibly discounted) cost functions (Bertsekas and Shreve, 1996). Starting from
the regression bounds in Munos and Szepesvari (2008), this work includes new results on the
error for the FVI approximation and a-priori performance guarantees. Additionally, novel
and tighter probabilistic error bounds for dynamic programming solutions of the reach-avoid
problem based on samples extraction are presented. As a comparison to the alternative techniques in the literature (Abate et al., 2010; Esmaeil Zadeh Soudjani and Abate, 2013), we
show the related techniques provide bounds that are valid deterministically, whereas the
proposed result yields tighter results in general that are valid with a certain (tunable) confidence. The outcomes lead to an approach providing controller synthesis with a certified
performance, which is relevant for safety-critical applications (Blom and Lygeros, 2006).
The proofs of the statements are included in the Appendix.

2. Probabilistic Reach-Avoid Problem over General Markov Processes

Definition 1 (General Markov process) A discrete-time general Markov process is comprised of:
• A continuous (uncountable) state space X ⊂ Rn ;
• An action space A = {a1 , . . . , am } consisting of a finite number of actions;
• A Borel-measurable stochastic kernel Tx , which assigns to each state-action pair x ∈ X
and a ∈ A a probability distribution Tx (· | x, a) over X .
We denote with (X , B(X ), P ) a probability structure on X , where
R B(X ) is the σ-algebra
associated to X and P is characterized as P (y ∈A|xR∈ X , a ∈ A) = A TxR(dy|x, a). We assume
that the stochastic kernels admit densities so that A Tx (dy | x, a) = A tx (y | x, a) dy.

Definition 2 (Markov policy) A Markov policy µ over horizon [0, Nt ] is a sequence µ =
(µ0 , µ1 , . . . , µNt −1 ) of universally measurable maps, µk : X → A, k = 0, 1, . . . , Nt − 1, from
the state space X to the action space A. The set of Markov policies is denoted as M.
3

Haesaert, Babuska and Abate

The evolution of the general Markov process is considered over a finite horizon k = 0, 1, . . . , Nt ,
with Nt ∈ N. Consider a discrete-time general Markov process, a Markov policy µ, a deterministic initial state x0 ∈ X and a finite time horizon Nt : an execution of the process
characterizes a state trajectory given as {xk |k = 0, 1, . . . , Nt }. The process evolves over
the product space (X )Nt +1 , which is again endowed with a (product) σ-algebra and allows
computing probability associated to events over trajectories – we denote this probability
by P, and further define the probabilities Px0 , Pµx0 as P conditioned on an initial state
and on an initial state and a policy, respectively. The state at the (k + 1)-st time instant, xk+1 , is obtained as a realization of the controlled Borel-measurable stochastic kernel
Tx (· | xk , µk (xk )). The model can be initialized according to an initial probability measure
P0 ∈ M (X ), where M (·) denotes the collection of probability distributions over a given set.
2.1 Probabilistic Reach-Avoid Problem: Definition
Let us define the probabilistic reach-avoid problem, also known as constrained reachability (Baier and Katoen, 2008), and provide its characterization. Consider a safe set
A ∈ B(X ), a target set K ∈ B(X ), and a finite time horizon Nt ∈ N. A given state trajectory {xk |k = 0, 1, . . . , Nt } verifies the reach-avoid property if it reaches the target set K
within the time horizon, while staying inside the safe set A. This property can be expressed
as
∃j ∈ [0, Nt ] : xj ∈ K ∧ ∀i ∈ [0, j − 1] : xi ∈ A \ K.
Let us now consider the probabilistic reach-avoid property for a general stochastic system,
defined as the probability that an execution associated with a fixed Markov policy µ ∈ M
and an initial condition x0 ∈ X reaches the target set K while avoiding X \ A. Formally,
rxµ0 (K, A) = Pµx0 {∃j ∈ [0, Nt ] : xj ∈ K ∧ ∀i ∈ [0, j − 1] : xi ∈ A \ K} ,

(1)

where the states x0 , x1 , . . . , xNt ∈ X are sampled via the stochastic kernel Tx under policy
µ. The formula contained in (1) can be written as a boolean expression using indicator
functions, which leads to an expectation over the state trajectories as




P
Qj−1
µ
µ
rx0 (K, A) = Ex0
i=0 1A\K (xi ) 1K (xj ) ,
j∈[0,Nt ]

where 1B (x) = 1 if x ∈ B, else it is equal to 0. The reach-avoid problem subsumes other
known problems widely studied in System and Control Theory and in Formal Verification,
such as that of reachability of set K, which is simply obtained by selecting A = X , or that
of invariance within a set B, which is characterized as the dual of the reachability problem
over set X \ B.
For a given policy µ, the time-dependent value function Wk : X → [0, 1], defined as
Wkµ (x) =



Eµ

P



j−1
Q



1A\K (xi ) 1K (xj ) xk = x ,

j∈[k+1,Nt ] i=k+1

4

Quantitative Approximations for Reach-Avoid over General Markov Processes

is the probability that the state trajectory {xk+1 , . . . , xNt }, starting from xk , will reach
the target set K within the time horizon [k, Nt ], while staying within the safe set A. This
function allows expressing the reach-avoid probability backward recursively, as follows.
Proposition 3 Given a policy µ = (µ0 , µ1 , . . . , µNt −1 ), define function Wkµ : X → [0, 1] by
backward recursion


µ
Wkµ (x) = Eµxk 1K (xk+1 ) + 1A\K (xk+1 )Wk+1
(xk+1 ) ,
and initialized with WNµt (x) = 0. Then for any initial state x0 ∈ X , the probabilistic reachavoid property rxµ0 (K, A) can be expressed as
rxµ0 (K, A) = 1K (x0 ) + 1A\K (x0 )W0µ (x0 ) .
Proof The proof follows (Summers and Lygeros, 2010, Lemma 4), where the above statement is proven for a value function Vkµ (x) = 1K (x) + 1A\K (x)Wkµ (x).

Notice that, while the probabilistic reach-avoid problem has been formulated above via DP
recursions, it hinges on a sum-multiplicative characterization which is non-standard: much
of the analytical and computational results in DP are formulated for additive (possibly
discounted) cost functions (Bertsekas and Shreve, 1996).
Rather than selecting and fixing a policy µ as done above, we now focus on the controller synthesis problem, which seeks the Markov policy µ∗ that maximizes the probabilistic
reach-avoid property, and which is such that rx∗0 (K, A) = supµ∈M rxµ0 (K, A). Let us emphasize that the optimization is over finite-action policies, which are however functions of
the continuous state space. The optimal policy can be characterized as follows.
Proposition 4 Define functions Wk∗ : X → [0, 1], by the backward recursions


∗
(xk+1 ) ,
Wk∗ (x) = max Eax 1K (xk+1 ) + 1A\K (xk+1 )Wk+1
a∈A

with xk+1 ∼ Tx (· | x, a) for k = Nt − 1, Nt − 2, . . . , 0, and initialized by WN∗ t (x) = 0. Then
for any initial state x0 ∈ X the optimal probabilistic reach-avoid property rx∗0 (K, A) can be
expressed as
rx∗0 (K, A) = 1K (x0 ) + 1A\K (x0 )W0∗ (x0 ).

Furthermore, µ∗k : X → A for k = Nt − 1, Nt − 2, . . . , 0, is such that ∀x ∈ X :


∗
µ∗k (x) = arg max Eax 1K (xk+1 )+1A\K (xk+1 )Wk+1
(xk+1 )
a∈A

and µ∗ = (µ∗0 , µ∗1 , . . . , µ∗Nt −1 ) is the optimal probabilistic reach-avoid Markov policy.

Proof See again (Summers and Lygeros, 2010, Theorem 6) and (Abate et al., 2008).

5

Haesaert, Babuska and Abate

For a given time horizon Nt , the computation of rx∗0 (K, A), as in Proposition 4, can be
seen as the application of Nt mappings. More precisely, let us define a dynamic programming
∗ , such that for all states x ∈ X , the function W ∗ : X → [0, 1] is
operator T as Wk∗ = TWk+1
k
defined as

∗
Wk∗ (x) = TWk+1
(x)


a
∗
= max Ex 1K (xk+1 )+1A\K (xk+1 )Wk+1
(xk+1 ) .

(2)

a∈A

The value of the optimal probabilistic reach-avoid property can be written as the composition of Nt mappings,

rx∗0 (K, A) = 1K (x0 ) + 1A\K (x0 ) TNt WN∗ t (x0 ).
2.2 Computation of the Reach-Avoid Probability

Notice that generally it is not possible to solve the above recursions exactly: in order
 to
∗
∗
determine the backwards iteration at a single point xi ∈ X , namely Wk (xi ) = TWk+1 (xi ),
one should exactly solve (2). The exact solution of (2) however is seldom analytical and
can possibly result in computationally expensive procedures. The absence of an analytical
representation for Wk∗ : X → [0, 1] leads to the use of approximation techniques, which can
be categorized in two families:
1. Numerical approximation techniques, which provide an approximation of the optimal probabilistic reach-avoid problem with actual error bounds. More precisely, for
a given error bound ∆ > 0 we seek a numerical scheme that obtains an approximation r̂x∗0 (K, A), which is such that |r̂x∗0 (K, A) − rx∗0 (K, A)| ≤ ∆. This approach
is taken in (Abate et al., 2010; Esmaeil Zadeh Soudjani and Abate, 2013), and the
scheme is prone to suffer from the curse of dimensionality, since it approximates a
general stochastic system with a Markov chain by partitioning the state space.
2. Probabilistic approximation techniques, which approximate the original problem with
probabilistic guarantees. The obtained approximation scheme r̂x∗0 (K, A) for rx∗0 (K, A)
depends on a finite number of samples or sampled paths of the underlying model. For a
given error bound ∆ > 0 and confidence 1−δ∆ , the probability that the approximation
is not close to the the optimal value can be bounded probabilistically as
P



r̂x∗0 (K, A) − rx∗0 (K, A) > ∆ ≤ δ∆ .

(3)

In this work, we newly pursue the second approach by focusing on results from the area
of learning, and in particular on algorithms for functional approximations. A learning approach is suitable for complex systems, such as general Markov processes, since it replaces
model-based evaluations by model-free, sample-based evaluations (Busoniu et al., 2010).
We adopt the Fitted Value Iteration scheme (FVI) (Munos and Szepesvari, 2008), a learning algorithm fit in particular for finite horizon settings.
In practice, bounds for probabilistic approximation methods (3) can be divided into two
groups. Firstly, general model-free and sample-free bounds, which provide a-priori guarantees on the achievable accuracy for a finite sample set. Though they can show convergence
in probability up to a bias term, their generality can render them conservative when used
6

Quantitative Approximations for Reach-Avoid over General Markov Processes

as a tool to assign an accuracy guarantee. Alternatively, model-based and sample-based
bounds: these bounds verify the accuracy of a dynamic programming scheme by drawing
samples of the model and using available information from the model, and from the specific
reach-avoid property under study. In the analysis of the algorithm these bounds can be perceived as complementary. A-priori and sample-free bounds are derived in Section 4 based
on model-free/distribution-free notions, whereas model-based and sample-based bounds are
given in Section 5.

3. Fitted Value Iteration
In this section we consider a learning algorithm that has been developed to solve additivecost optimal control problems, and adapt it to the reach-avoid optimal control setting.
Known as FVI Munos and Szepesvari (2008), the algorithm extracts a finite number of
samples from the underlying model to numerically approximate the value recursions in
(2). More precisely, the scheme generalizes the information gathered from the samples to
approximate the “exact” optimal value function Wk∗ as Ŵk∗ in two steps: first by estimating
Wk∗ over a finite number of states, thereafter fitting the analytical function Ŵk∗ ∈ W to the
estimate. We define W to be a strict subset of B(X ; 1), the class of measurable functions
defined over X , lower bounded by 0 and upper bounded by 1.
Algorithm 1: Sample generation for the FVI algorithm
Given a safe set A, a target set K, a time horizon Nt , and distribution η, generate N base
points, and M samples at each base point as follows:
1. Draw N base points xik 1≤i≤N from the distribution η in A \ K;
2. Draw M samples at each base point xik andat each
 actions a ∈ A from the stochastic
kernel Tx , and denote the set of samples as xi,a,j
k+1

1≤j≤M

.

The FVI algorithm employs samples that are generated by the underlying model. Samples referred to as “base points” are taken from a chosen distribution η, and additional
samples are drawn at the base points from the transition kernel. Algorithm 1 summarizes
the sample generation. Let us remark that at each iteration k = Nt − 1, . . . , 0, a new set
of samples is generated and used. We denote as “sample complexity” the cardinality of the
sample generation, namely N and M .
At each (backward) iteration k = Nt − 1, . . . , 0, the algorithm executes two steps in
order to approximate the exact recursion in (2):
∗ )(xi )
1. The first step consists of estimating the value of the backward mapping (TŴk+1
k
∗ )(xi ) is estimated by an empirical operator
at N base points xik . The recursion (TŴk+1
k
T̂ as follows:
M
 i
1 X
i,a,j
i,a,j
∗
∗
1K (xi,a,j
(4)
T̂Ŵk+1
(xk ) = max
k+1 ) + 1A\K (xk+1 )Ŵk+1 (xk+1 ).
a∈A M
j=1

xi,a,j
k+1

Here
represent M independent and identically distributed realizations obtained

from Tx · | xik , a . Hence, j ∈ [1, M ], i ∈ [1, N ], and a ∈ A. For an increasing number
7

Haesaert, Babuska and Abate

∗
i
∗
i
of samples xi,a,j
k+1 , the estimate (T̂Ŵk+1 )(xk ) converges to (TŴk+1 )(xk ) with probability
1, by the law of large numbers.

2. In the second step, function Ŵk∗ ∈ W is estimated as the solution of
Ŵk∗ = arg min

w∈W

N
X
i=1

p

∗
(xik ) .
w(xik ) − T̂Ŵk+1

(5)

We assume that the argument of the minimum belongs to the function class (this fact
will be further discussed below). The base points (xik )1≤i≤N are independently drawn
from a distribution η supported over the set A \ K. The power factor p ≥ 1 is a given
positive number. Given
a function w(x) and an increasing value of N , the summation
R
∗ (x)|p η(x)dx, by the law of large numbers. This
in (5) converges to A\K |w(x)− T̂Ŵk+1

leads (cf. Section 4) to the convergence of the argument Ŵk∗ to the optimal fit that
∗
minimizes the distance between T̂Ŵk+1
and the functions w ∈ W with respect to the
p-norm, weighted by a distribution with density η and supported over the set A \ K,
namely
∗
kŴk∗ − T̂Ŵk+1
kp,η

=

Z

|Ŵk∗ (x)
A\K

p1

∗
− T̂Ŵk+1
(x)|p η(x)dx

.

The overall FVI algorithm is summarized in Algorithm 2. The iterations are initialized
as ŴN∗ t (x) = 0, x ∈ X , and updated over the functions ŴN∗ t −1 , ŴN∗ t −2 , . . . , Ŵ1∗ . Finally the
value function at k = 0 is approximated at the initial condition x0 with a sample-based
integration, similar to (4), using M0 independent and identically distributed realizations of
Tx (· | x0 , a) for all a ∈ A.
Algorithm 2: Fitted Value Iteration algorithm
Given an initial condition x0 ∈ X , a safe set A, a target set K, a time horizon Nt , a set of
N base points and of M samples at each base point (for each iteration k), a number p and a
distribution η, perform:
1. Initialize ŴN∗ t (x) = 0, ∀x ∈ X ;
2. For k = Nt − 1 to 1 do
(a) Collect samples (cfr. Algorithm 1);
 i
∗
∗
(b) Estimate TŴk+1
as T̂Ŵk+1
(xk ) (cfr. (4))
(c) Find the function that minimizes the empirical p-norm as
Ŵk∗ = arg min

w∈W

N
X
i=1

p

∗
(xik ) ;
w(xik ) − T̂Ŵk+1

3. Collect M0 samples for the single initial condition x0 and for every action a ∈ A, and
estimate Ŵ0∗ (x0 ) as in step (2)(b);
4. Return reach-avoid probability
r̂x∗0 (K, A) = 1K (x0 ) + 1A\K (x0 )Ŵ0∗ (x0 ), ∀x0 ∈ X .

8

Quantitative Approximations for Reach-Avoid over General Markov Processes

Remark 5 (Approximately Optimal Policy) The FVI algorithm can be extended to
include the synthesis of a policy a = µ̂k (x). At every iteration in Algorithm 2, first the policy
is estimated at all the base points xik , as the argument of (2)(b). Secondly a classification
algorithm is used, providing an approximately optimal policy µ̂k : X → A for each k.

4. A-Priori Probabilistic Error Bounds
Let us recall the accuracy of the FVI algorithm as in (3): we say that the FVI algorithm
has an accuracy ∆, with a confidence 1 − δ∆ , if the probability that the error made by the
approximate solution r̂x∗0 (K, A) is larger than ∆, is upper-bounded by δ∆ . We explicitly
quantify the accuracy in (3) and analyse it in two steps: first by computing a bound on
the error of a single iteration (Sec. 4.1); then by studying the propagation of the singlestep error over multiple iterations (Sec. 4.2). Notice that this section provides a-priori
bounds which are model- and distribution-free, hence computable before applying the FVI
algorithm. Alternatively, a-posteriori error bounds based on an additional sample set are
proposed in Section 5.
4.1 Error Bounds on a Single Iteration of the FVI
∗
: X → [0, 1] be an unknown map, and consider a function class W ⊂ B(X ; 1).
Let TŴk+1
Recall that at each iteration k = Nt − 1, Nt − 2, ..., 1, the objective of the learning algorithm
∗
with respect to the following weighted,
is to find a function w ∈ W that is close to TŴk+1
p-norm:
1
Z
p
p
∗
∗
(6)
w − TŴk+1 η(x)dx .
kw − TŴk+1 kp,η =
X

∗ k
∗
6∈ W, the optimal approximation inf w∈W kw − TŴk+1
Notice that if TŴk+1
p,η provides
only a lower bound on this error: the presence of this non-zero bias error indicates that
the FVI scheme is not asymptotically consistent, namely that the error does not converge
to zero for an increasing sample size. Of interest to this work, a general upper bound on
∗ k
inf w∈W kw − TŴk+1
p,η is derived as

dp,η (TW, W) = sup inf kf − Tgkp,η .
g∈W f ∈W

(7)

In (Munos and Szepesvari, 2008) the bias dp,η (TW, W) is referred to as the inherent Bellman
error of the function space W.
As discussed, the FVI algorithm employs empirical estimates, given in (4)-(5), of the
quantity in (6). Therefore, the single step error hinges both on the inherent Bellman error,
∗
and on the deviations caused by using estimates of the recursion step TŴk+1
over the base
points (cfr. Section 4.1.1) and of the norm k · kp,η as the integral in (6) (cfr. Section 4.1.2).
The error contributions depend on the number of samples used (N, M ) and on the capacity
of the function class W (Section 4.1.2, and Appendix B), whereas they do not depend on the
distribution η, nor on the stochastic state transitions characterizing the model dynamics:
as such the bounds are general and “distribution-free”.
In the following subsections, two lemmas are derived, which are necessary to obtain a
general upper bound for the single step error. First (Section 4.1.1), the error introduced by
9

Haesaert, Babuska and Abate

∗
using the estimation T̂Ŵk+1
of the recursion step is bounded using Hoeffding’s inequality
(Hoeffding, 1963). Then in Section 4.1.2 the maximal deviation of the empirical evaluation
of the integral in (6) is bounded using methods from Statistical Learning Theory (Vapnik,
1998).
∗
4.1.1 Accuracy of the Estimation of TŴk+1
 i
 i
∗
∗
Recall that the estimate T̂Ŵk+1
(xk ) of the exact recursion TŴk+1
(xk ) uses, for a given
i
state-action pair (xk , a), the individual Monte-Carlo estimates
M
1 X
i,a,j
i,a,j
∗
1K (xi,a,j
k+1 ) + 1A\K (xk+1 )Ŵk+1 (xk+1 )
M
j=1

i
h
∗ (x
of Eaxi 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
k+1 ) . Since the cardinality of the action space A
k
is finite, a bound on the error of the estimates for a given state-action pair can lead to a
bound on the error over the states xik : we elaborate on this idea next.
i,a,j
i,a,j
∗
The M random quantities for 1 ≤ j ≤ M , 1K (xi,a,j
k+1 ) + 1A\K (xk+1 )Ŵk+1 (xk+1 ), are
obtained via independent and identically distributed realizations over the closed interval
[0, 1]. Hoeffding’s inequality (Hoeffding, 1963) leads to an upper bound on the deviation of
the estimate from the expected value as follows:
n
h
i
∗
P Eaxi 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
(xk+1 )
k

M
i
o
1 Xh
2
i,a,j
i,a,j
∗
(x
)
≥ 1 − 2e−2M (ǫ1 ) ,
1K (xi,a,j
)
+
1
(x
)
Ŵ
≤
ǫ
−
1
A\K
k+1 k+1
k+1
k+1
M
j=1

where ǫ1 is the bound on the error.
We can then provide
 i a lower bound on the probability that the deviation incurred by the
∗
quantity T̂Ŵk+1
(xk ) is bounded by ǫ1 via the joint probability of |A| independent events,
as follows:
o
n
∗
∗
P TŴk+1
(xik ) − T̂Ŵk+1
(xik ) ≤ ǫ1
h
i
Y n
∗
(xk+1 )
≥
P Eaxi 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
a∈A

k

M
o
1 X
i,a,j
i,a,j 
∗
−
.
1K (xi,a,j
)+1
(x
)
Ŵ
(x
)
≤
ǫ
1
A\K
k+1 k+1
k+1
k+1
M
j=1

Let us now extend the above probabilistic bound for the error computed at a single point
xik to a bound for the error over all the N base points: we can express this bound via an
empirical p-norm defined over the base points, as follows:
∗
∗
kp,η̂ =
− T̂Ŵk+1
kTŴk+1

N
p 1/p
1 X
∗
∗
.
(xik )
(xik ) − T̂Ŵk+1
TŴk+1
N
i=1

This leads to the following result.
10

(8)

Quantitative Approximations for Reach-Avoid over General Markov Processes

Lemma 6 For a given error bound ǫ1 and sample complexity N and M , the estimation
error can be probabilistically bounded as follows:
o
n
∗
∗
(9)
kp,η̂ ≤ ǫ1 ≥ 1 − δ1 ,
− T̂Ŵk+1
P kTŴk+1
2

2

where δ1 = 1 − (1 − 2e−2M (ǫ1 ) )|A|N , as long as 0 < 2e−2M (ǫ1 ) ≤ 1.

Notice that for an increasing number of samples M , by Lemma 6 the empirical norm of the
error as in (8) is less than ǫ1 with a probability that increases to 1.
4.1.2 Accuracy of the Empirical Norm
∗
: X → [0, 1] be an unknown function and η a probability measure on X ,
Let TŴk+1
∗
with respect
η ∈ M (X ). The objective is to find a function w ∈ W that is close to TŴk+1
to the following expected loss:
R
∗
∗ (x)|p η(x)dx.
inf kw−TŴk+1
kpp,η = inf |w(x)−TŴk+1
w∈W

w∈W X

This section provides a bound for the error originating from the use of a finite number of
samples to evaluate the loss: this empirical loss is defined as
kw −

∗
kpp,η̂
TŴk+1

N
1 X
∗
=
(xik )|p ,
|w(xik ) − TŴk+1
N
i=1

for a given set of N random variables drawn independently over A \ K according to xik ∼ η.
Let us express a probabilistic bound on this error that holds uniformly over all functions
w ∈ W as follows:

∗
∗
P sup kw − TŴk+1
kpp,η−kw−TŴk+1
kpp,η̂ ≥ ǫp2 ≤ δ2 .
w∈W

Observe that since the expected and the empirical losses can be reformulated respectively
as the mean and the empirical mean of a loss function, defined informally as f (x) =
∗ (x)|p , the above problem can be framed as the uniform convergence of a
|w(x) − TŴk+1
standard learning problem (Haussler, 1992; Pollard, 1984; Vapnik, 1998). Resorting to related literature, results on bounds for the error probability of the regression of real-valued
functions employ capacity concepts (including Rademacher averages, covering numbers, and
pseudo dimensions) of a function class (Bartlett et al., 2005; Hoeffding, 1963). We focus on
pseudo dimensions to deal with the capacity (or the complexity) of the function class. By
results in (Haussler, 1992) and (Pollard, 1984), we obtain the following uniform convergence
bound.

Lemma 7 Let W be a set of finitely parameterized, measurable functions on X taking values
in the interval [0, 1] with pseudo dimension dimp (W) = d < ∞. Let (xik )1≤i≤N be generated
by N independent draws according to any distribution η on A \ K and let p ≥ 1. Then for
any ǫ2 > 0 we have that




2p
32e d − Nǫ2
∗
∗
128 .
e
(10)
P sup kw − TŴk+1
kpp,η − kw − TŴk+1
kpp,η̂ ≥ ǫp2 ≤ 4e(d + 1)
ǫp2
w∈W
The characterization and computability of the pseudo dimension of a function class is given
in Appendix B.
11

Haesaert, Babuska and Abate

4.1.3 Global Accuracy of Single Iterations
The error bounds introduced in Lemmas 6 and 7, together with the inherent Bellman error,
yield an upper bound on the error introduced at each iteration, which is recapitulated
in the following statement. Since the bounds are independent of the distributions η and
Tx (· | x, a), they are in fact distribution-free bounds (Bartlett et al., 2005).
Theorem 8 Consider a reach-avoid property defined over a general stochastic system with
continuous state space X and finite action space A. Let A ⊂ X and K ⊂ X be Borel
measurable sets and fix p ≥ 1, the distribution η ∈ M (A \ K) and W ⊂ B(X ; 1). Pick any
∗
∗
Ŵk+1
∈ B(X ; 1) and let T̂Ŵk+1
and Ŵk∗ be calculated using (4) and (5). For given upper
bounds on

∗
∗
P kTŴk+1
− T̂Ŵk+1
kp,η̂ > ǫ1 ≤ δ1 , and

∗
∗
P sup kw − TŴk+1
kpp,η − kw − TŴk+1
kpp,η̂ > ǫp2 ≤ δ2 ,

(11)
(12)

w∈W

the bound on a single step update is as follows:
n
o
∗
∗
P kŴk∗ − TŴk+1
kp,η > dp,η (TŴk+1
, W) + ǫ ≤ δ1 + δ2 ,

(13)

∗ , W) =
where the optimal approximation is given as a biasing term defined as dp,η (TŴk+1
∗ k , and the error ǫ is given as ǫ = 2ǫ + 2ǫ .
inf w∈W kw − TŴk+1
p,η
1
2

Note that the statement assumes that Ŵk∗ is the unique solution to the optimization problem
in (5). This strict assumption can be weakened by adding a tolerance term to the theorem.
∗ , W) admits the inherent Bellman error d (TW, W) in (7) as a
The quantity dp,η (TŴk+1
p,η
general upper bound.
4.2 Error Propagation and Global Error Bounds
We express a global bound on the accuracy of the FVI algorithm r̂x∗0 (K, A) − rx∗0 (K, A)
by using the probabilistic bounds (derived
 on the errors introduced by the
 in Section 4.1)
approximate one-step mappings | T̂Ŵ1∗ (x0 ) − TŴ1∗ (x0 )|, as well as kŴ1∗ − TŴ2∗ kp,η ,
kŴ2∗ − TŴ3∗ kp,η , . . . , kŴN∗ t −1 − TŴN∗ t kp,η , and by propagating the error introduced by each
single iteration to the successive value iterations, as done in the next statement. More
precisely, in the following lemma we show that the deviation of the approximate value
function Ŵk∗ from the optimal value function TNt −k WN∗ t can be expressed as a function
of this deviation at step k + 1, plus the approximation error introduced at the (Nt − k)th
iteration (which has been bounded above). Recall that the optimal value functions can be
written as Wk∗ = TNt −k WN∗ t and that by definition WN∗ t = ŴN∗ t .
Lemma 9 Let η be the density of a probability distribution with support on A \ K and tx
be the density function of the stochastic kernel Tx . Then
1

∗
∗
− TNt −(k+1) ŴN∗ t
kp,η + B p Ŵk+1
kŴk∗ − TNt −k ŴN∗ t kp,η ≤ kŴk∗ − TŴk+1

12

p,η

,

Quantitative Approximations for Reach-Avoid over General Markov Processes

for k = 0, 1, . . . , Nt − 1, and where B is defined as
Z max tx (xk+1 | xk , a) η(xk )
a∈A
dxk .
B = sup
η(xk+1 )
xk+1 ∈A\K A\K
Putting all the pieces together, the following theorem provides an expression for the global
FVI error bound as the accumulation of the errors from the single iterations over the whole
time horizon.
Theorem 10 Consider a reach-avoid problem defined on a Markov process with a continuous state space X and a finite action space A. The optimal reach-avoid probability rx∗0 (K, A)
for a given target set K, safe set A, initial state x0 , and time horizon Nt , is approximated
by the quantity r̂x∗0 (K, A) obtained with the FVI Algorithm in Algorithm 2, which has an
accuracy of ∆ and a confidence δ∆ , as stated in (3), if the following holds:
 NX
t −1
k−1
∗
P B0
kp,η
B p kŴk∗ − TŴk+1

(14)

k=1

+

|Ŵ0∗ (x0 )

−

TŴ1∗ (x0 )|

where B is given in (14), B0 is defined as
tx (x1 | x0 , a)
,
B0 = sup max
η(x1 )
x1 ∈A\K a∈A


> ∆ ≤ δ∆ ,

(15)

(16)

and where η is the density of a probability distribution supported on A \ K and tx is the
density of the stochastic kernel Tx of the given Markov process.
Remark 11 Notice that the scaling factor B has a significant influence on the error propagation. It is related to the notion of concentrability of the future-state distribution (Munos and Szepesvari,
2008; Farahmand et al., 2010). If B > 1 the error of the algorithm will increase exponentially with the time horizon, whereas if B < 1 the accuracy of the algorithm will depend
mostly on the errors in the last few iterations (backwards in time). It B expresses the maximal concentration of the dynamics (relative to η) over the relevant set A \ K after one
transition starting from distribution η.
The case study discussed in Section 6 displays a choice of a density function η leading
to bounded values for B and B0 , respectively. The relation between the scaling factor B
and the model dynamics for the considered Markov process is also analyzed.
4.2.1 Discussion on the Global Error Bounds
Suppose that we require that the error in the estimation of rx∗0 (K, A), with a confidence at
least equal to α, is less than ∆, as per (3). Let us assume that there exist positive values
ǫ0 , ǫ1 , ǫ2 such that the approximation error ∆ can be split up into individual bounds based
on Theorem 8 and 10 as follows:
∆ =B0

NP
t −1
k=1

B

k−1
p

dp,η (TW, W) + 2B0
13

NP
t −1
k=1

B

k−1
p

(ǫ1 )

Haesaert, Babuska and Abate

+ 2B0

NP
t −1

B

k−1
p

(ǫ2 ) + ǫ0 ,

k=1

accounting for, respectively, the inherent Bellman error (Section 4.1); the error on the
∗
estimation of TŴk+1
(Section 4.1.1); the error on the empirical norm (Section 4.1.2); and
the error related to |Ŵ0∗ (x0 ) − TŴ1∗ (x0 )|. We compute a lower bound on the confidence
that the approximation error is bounded by ∆ by upper bounding the complementary event.
Using a union bounding argument on the probability of invalidating any of the bounding
terms in the above equation, we obtain that the probability that the overall error is larger
than ∆ is upper bounded with δ∆ , as

2
(17)
δ∆ =0 + (Nt − 1) 1 − (1 − 2e−2M (ǫ1 ) )|A|N


 d Nǫ2p
2
e− 128
+ (Nt − 1) 4e(d + 1) 32e
ǫp
2

+ 1 − (1 −

2
2e−2M0 (ǫ0 ) )|A| ,

with respectively the inherent Bellman error (Section 4.1), the confidence terms from Lemma
6 with sample complexity M and N , from Lemma 7 with dimp (W) = d < ∞, and from
2
Lemma 6 with M = M0 and N = 1. Eqn. (17) holds as long as 0 < 2e−2M (ǫ1 ) ≤ 1 and
2
0 < 2e−2M0 (ǫ0 ) ≤ 1. We can observe that it is possible to find finite values for N ,M ,M0
such that 1 − α > δ∆ . Moreover for each choice of positive ǫ0 , ǫ1 , ǫ2 > 0 and confidence
0 ≤ α < 1, the necessary number of samples can be upper bounded by polynomials in
1 1 1
1
ǫ0 , ǫ1 , ǫ2 and 1−α , as follows:
l
 2p
N = 128 ln(4e(d + 1)) + d ln(32e) ǫ12
 2p  m
 2p  
ln ǫ12 + 128 ǫ12
ln δ12 ,
+ 128dp ǫ12
m
l  2 
ln(2|A|) + ln( δ11 ) + ln(N ) ,
M = 21 ǫ11
l  2 

m
ln 2|A| + ln δ10
,
M0 = 21 ǫ10

(18a)

(18b)
(18c)

with positive parameters δ0 , δ1 , δ2 > 0 such that 1 − α = δ0 + (Nt − 1)δ1 + (Nt − 1)δ2
(derivation in Appendix F).
Note that the above accuracy does not depend on the dimensionality n of the state space.
Therefore the accuracy for models with higher state space dimension will directly depend on
the complexity of the function class employed to approximate given value functions. This
is unlike standard grid-based numerical approximation techniques, such as that proposed
in (Esmaeil Zadeh Soudjani and Abate, 2013), which are known to break down over models
with large state-space dimensionality.
Let us add a few comments on the dependency of the accuracy from several design variables
of the FVI algorithm. Firstly, the choice of function class affects both the inherent Bellman
error and the pseudo dimension: while the former gives a measure of how well the the
function class W can represent the value functions Wk∗ , the latter is directly related to
the complexity of the function class. The objective is to obtain a low complexity function
14

Quantitative Approximations for Reach-Avoid over General Markov Processes

class that is capable to accurately fit the given value functions. A good accuracy can be
hard to attain when a bad choice of the function class leads to both a large bias (due
to the inherent Bellman error) and to a large number of samples. Secondly, the sample
distribution η defines, together with the state transitions, the scaling factors B and B0 . In
order to minimize the error propagation caused by B, the distribution η should be “aligned”
with the model dynamics characterized by the density of the transition kernel. Finally, the
parameters N, M, M0 follow from the required accuracy demands, which are reformulated
as polynomial functions depending on ǫ0 , ǫ1 , ǫ2 and δ0 , δ1 , δ2 , e.g. as in (18a).
With regards to the single-step errors, Lemmas 6 and 7 determine a bound uniformly over
the whole function class and for any possible probability distribution. The used distributionfree notions lead to conservative bounds (Bartlett et al., 2005), which then result in a large
set of required samples. Notice however that the construction allows to compute the bounds
a-priori, before any sample is drawn from the system.
In conclusion, the formal probabilistic bounds on the error made by the approximation
algorithm show that the algorithm converges in probability to the best approximation for
an increasing cardinality of the samples.

5. Sample-Based Error Bounds
In this section, a probabilistic bound on the error of the approximated reach-avoid probability is developed according to a model- and sample-based philosophy. This bound can be
computed after the reach-avoid probability has been obtained via dynamic programming
as time-dependent, approximate value functions Ŵk∗ for 0 ≤ k ≤ Nt − 1. The obtained
bounds are not only sample dependent but also distribution dependent, since knowledge of
the transition kernel is necessary to compute scaling factors such as (14). Throughout the
section it is assumed that the used samples are not correlated with Ŵk∗ , in other words if
the estimated value functions Ŵk∗ are a result of a sampled-based optimization, then the
samples used for the bounds in this section are drawn anew and independently.
The probabilistic bound on the accuracy P { | r̂x∗0 (K, A) − rx∗0 (K, A) | > ∆ } ≤ δ∆ , as
given in (3), and computed now in a sample-based manner, includes an empirical estimate of
the quantity | r̂x∗0 (K, A) − rx∗0 (K, A) |. This sample-based estimate is computed as follows:
a. collect samples (xi )1≤i≤Ñ according to (1) in Algorithm1 (with N = Ñ ), and subsequently use (2) with M = M̃ to draw both (y1i,a,j )1≤j≤M̃ and (y2i,a,j )1≤j≤M̃ ;
b. estimate the single step error (19) for each k;
c. estimate the bias (21) for each k;
d. compute the multi-step error as a propagation and a composition of the estimates in
(b.) and (c.).
The single step error is estimated as
∗
Ŵk∗ − T̂Ŵk+1

=
1,η̃

1
Ñ

Ñ
P

i=1

∗ (xi )
Ŵk∗ (xi ) − max T̂a1 Ŵk+1
a∈A

(19)

for all 1 ≤ k ≤ Nt −1. The term on the right is the empirical 1-norm and can be written as a
1-norm with weighting η̃, which is the empirical distribution of η resulting from (xi )1≤i≤Ñ .
15

Haesaert, Babuska and Abate

Let the operator T̂aα be, for α = 1, 2,
∗ (xi ) = 1
T̂αa Ŵk+1
M̃

M̃
P

∗ (y i,a,j ).
1K (yαi,a,j )+1A\K (yαi,a,j )Ŵk+1
α

(20)

j=1

The estimate in (19) is biased due to the maximization over the action space, therefore
as a second step we estimate a bound on this bias. The combination of the two sample
sets (y1i,a,j )1≤j≤M̃ and (y2i,a,j )1≤j≤M̃ , for each xi and a, allows us to estimate this bias for
1 ≤ k ≤ Nt − 1 as
∗
∗
− T̂a2 Ŵk+1
max T̂a1 Ŵk+1
a∈A

1,η̃

.

(21)

In the following theorem, an expression for the bound on P { | r̂x∗0 (K, A) − rx∗0 (K, A) | >
∆ } ≤ δ∆ is derived, employing the error propagation technique first used in Section 4.2, the
estimates of the single step error above, and the bias (21) in combination with Hoeffding’s
inequality (Hoeffding, 1963).
Theorem 12 Consider a reach-avoid problem defined on a Markov process with a continuous state space X and a finite action space A. The optimal reach-avoid probability rx∗0 (K, A)
for a given target set K, safe set A, initial state x0 , and time horizon Nt , is approximated
by the quantity r̂x∗0 (K, A) obtained with the FVI Algorithm in Algorithm 2, which has an
accuracy of ∆ and a confidence δ∆ , as stated in (3), if the following holds:


NP
t −1
k−1
∗
∗
a
∗
a
∗
∆ = B0
B
Ŵk − T̂Ŵk+1
+ maxa∈A T̂1 Ŵk+1 − T̂2 Ŵk+1
+ B0 ǫ + ǫ 0 ,
1,η̃

1,η̃

k=1

(22a)

2
−2 Ñǫ2
L

δ∆ = e

− δ0 ,

(22b)

PNt −1 k−1
with L = 2 k=1
B
, and sample sizes M̃ and Ñ according to the sample sets drawn
according to the distribution η. Equation (22a) includes the estimated error as a combination
of (19) and (21). The scaling factors B and B0 are computed as in (14) and (16) for the
same sampling distribution η. The factors δ0 and ǫ0 are computed as in Lemma 6 with
M = M0 and N = 1.
The accuracy ∆ depends on two terms, the propagation of the estimated single-step and
bias errors over the time horizon up to k = 1, and the estimation errors B0 ǫ + ǫ0 for k = 0,
related to the confidence δ∆ .
Suppose that a close-to-optimal policy is given, for example a policy as detailed in
Remark 5 and computed from the series of estimated value functions Ŵk∗ . Then we know
∗
∗
that rx∗0 (K, A) ≥ rxµ̂0 (K, A), therefore a lower bound on the value of rxµ̂0 (K, A) is also a
lower bound on rx∗0 (K, A). Note that for a policy µ, the closed-loop Markov process is time
dependent. This allows us to estimate rxµ0 (K, A) directly from traces of this autonomous
Markov process. The deviation of this empirical mean can be bounded probabilistically
using Hoeffding’s inequality. Additionally an upper bound on the deviation |r̂x∗0 (K, A) −
rxµ0 (K, A)| can be computed. The combination of the bound in Theorem 12 and of the bound
on |r̂x∗0 (K, A) − rxµ0 (K, A)| provide a bound on the performance deviation of rxµ0 (K, A): the
16

Quantitative Approximations for Reach-Avoid over General Markov Processes

triangle inequality leads to |rx∗0 (K, A) − rxµ0 (K, A)| ≤ |rx∗0 (K, A) − r̂x∗0 (K, A)| + |r̂x∗0 (K, A) −
rxµ0 (K, A)|.
In comparison to the a-priori bound derived in Section 4, the sample-based bounds do
not depend on the inherent Bellman error and can be shown to be less conservative in
general. Moreover, they provide insight into the accuracy of the iterations steps. However,
they give no information about the expected convergence of the algorithm, and they can
only be computed after a run of the algorithm. Similarly to the a-priori bounds, they do
not depend on the dimensionality of the state space and are expected to scale better than
those used for grid-based approaches such as (Esmaeil Zadeh Soudjani and Abate, 2013).

6. Case Study and Numerical Experiments
We consider a case study from the literature (Fehnker and Ivančić, 2004), where the goal
is to maximize the probability that the temperature of two interconnected rooms, while
staying within a comfortable range, reaches a smaller target range within a given finite time
horizon. The temperature can be affected using local heaters. A reach-avoid problem is set
up by selecting as the safe set A = [17.5 22]2 , as the target set K = [19.25 20.25]2 , and
a fixed time horizon Nt = 10. The case study was implemented in Matlab R2013b on a
notebook with 2.6 GHz Intel Core i5 and 16 GB of RAM.
6.1 Model
The dynamics of the temperature in the two rooms is described by a Markov model, with
the temperature of the rooms making up the state space X = R2 , and where the possible
configurations {OF F, ON } = {0, 1} of the two heaters form the finite action space A. Hence
A = {0, 1} × {0, 1}, and as an example the action related to the first heater in the ON mode
and the second in the OFF one is given as a = [1 0]T ∈ A. The dynamics at discrete time
k is characterized by the following stochastic difference equation:
xk+1 = Axk + Ba + C + nk , where
i
h
i
h
 c1
1−b1 −a1,2
a1,2
b1 xa
A=
,
C
=
b2 xa , and B = 0
a2,1
1−b2 −a2,1

(23)
0
c2



,

and with the following parameters: xa is the ambient temperature (assumed to be constant),
bi ≥ 0 is a constant for the average heat loss rate of room i to the environment; aij ≥ 0 is
a constant for the average heat exchange rate of room i to room j 6= i; ci ≥ 0 is a constant
for the rate of heat supplied by the heater in room i. The parameters re instantiated
as b1 = 0.0375, c1 = 0.65, xa = 6, b2 = 0.025, c2 = 0.6, and aij = 0.0625. The noise
process nk is a realization of zero-mean Gaussian random variables with covariance ν 2 I2×2
(2-dimensional identity matrix I2×2 ) and ν = 0.5. Let N (· | µ, Σ) be a 2-dimensional
multivariate normal distribution over (X , B(X )) with mean µ and covariance matrix Σ,
then the stochastic kernel Tx is given as
Tx (· | x, a) = N · | Ax + Ba + C, ν 2 I2×2
17



(24)

Haesaert, Babuska and Abate

and characterises the probability distribution of the stochastic transitions in (23). The
stochastic kernel (24) admits the probability density
T −1
1
1
e(− 2 (y−µ̄) Σ (y−µ̄)) ,
tx (y | x, a) = p
|Σ|(2π)2

(25)

where | · | denotes the determinant of a matrix, and as before the covariance matrix equals
Σ = ν 2 I2×2 and the mean value is µ̄ = Ax + Ba + C.
6.2 Application of the Fitted Value Iteration Algorithm
The FVI scheme is implemented as in Algorithm 2, and approximates the solution of the
∗ (x ) by T̂Ŵ ∗ (x ),
reach-avoid problem. We obtain an approximation of rx∗0 (K, A) = T10 Ŵ10
0
0
1
∗
∗
∗
∗
while using the auxiliary functions Ŵ9 , . . . , Ŵ1 to approximate TŴ10 , . . . , TŴ2 in the FVI
scheme – equivalently, function Ŵk∗ approximates TNt −k ŴN∗ t . For a given temperature xk at
time instant k, the function Ŵk∗ (xk ) gives the approximate probability that the consecutive
temperature values xk+1 , . . . , xNt will reach the temperature range [19.25, 20.25]2 within
Nt − k time steps, while staying inside the safe set [17.5, 22]2 .

In order to apply the FVI algorithm, we select a uniform distribution η over A \ K to
sample from, then select a function class W and a value for p ≥ 1 to solve (5). We consider
a function class W composed of Gaussian radial basis function (RBF) neural networks
with 50 RBFs with a uniform width of 0.7. The neural network toolbox of Matlab is
used to solve the regression problem in (5) as a least-square problem (with p = 2). A
neural network with a single layer of hidden units of Gaussian type radial basis functions
is proved to be a universal approximator for real-valued functions (Hartman et al., 1990).
Furthermore the pseudo dimension of an artificial neural network with W free parameters
and k hidden nodes has been upper bounded by O(W 2 k2 ) (Karpinski and Macintyre, 1997;
Anthony and Bartlett, 1999). This means that for any desired precision the required number
of samples is bounded by a polynomial in the number of hidden nodes.

The following quantities are obtained for the sample complexities: N = 600, M = 103 ,
M0 = 103 . The approximate value functions for Ŵ9∗ , Ŵ5∗ , and Ŵ1∗ are displayed in Fig.
1. On the top plots, a point on the state space is associated with a probability for the
reach-avoid property over the given time horizon. At the bottom, the contour plots (level
sets) characterize the set of points that verify the reach-avoid property with a probability
at least equal to the given level.
Fig. 2 displays a suboptimal policy µ̂∗ that is obtained via the FVI algorithm as discussed in Remark 5, by employing the tree classification method ClassificationTree.fit
of Matlab. Observe that policy µ̂∗9 for k = 9 is not accurate over the flat regions of Ŵ9∗
(corresponding to the blue spots in the left side of Fig. 2 - left plot), which are far away
from the reach set K. Since the average heat loss rate of room 1 is the highest, we expect
that the heating should be turned on relatively longer. Fig. 2 confirms this, i.e. the red
(ON,ON) region is not square-shaped as the heaters stay ON for higher temperatures in
room 1 than in room 2.
18

Quantitative Approximations for Reach-Avoid over General Markov Processes

Ŵ5∗ (x)

Ŵ1∗ (x)

1

1

1

0.5

0.5

0.5

Ŵ9∗ (x)

18

19

0.6
20
x1

x2

x2

0.8

22

Ŵ1∗ (x)
0.8

21
20

22

18
22

0.8
0.6
0.4

20
x1

8

20
18 x1

0.8
0.6

18

22

0.

19
18

20
x2 18

0.6

18

3

0.

0
22

Ŵ5∗ (x)

0.4
0.2

0.1

20

20
18 x1

0.8
0.6
0.4

0.4

19

21

0.3

0.1

0.3
0.2
0.1

20

2
10.
0. 0.4

21

22

0.2

22

20
x2 18

22

x2

20
18 x1

0
22

0.8

20
x2 18

22

0.4
0.6

0
22

0.6

Ŵ9∗ (x)

18

20
x1

22

Figure 1: Function approximations Ŵ9∗ , Ŵ5∗ , and Ŵ1∗ – level sets (top) and contour plots
(bottom). The function approximations of the value functions are obtained with the FVI
algorithm as in Algorithm 2, using a radial basis function (RBF) neural network with 50
radial basis functions and a given width of 0.7, a uniform sampling distribution η over A\K
and M = 103 and N = 600. The safe set A is [17.5 22]2 , and reach set K = [19.25 20.25]2 .
∗ , T5 W ∗ , and T9 W ∗ over the set A \ K. In the
The Ŵ9∗ , Ŵ5∗ , and Ŵ1∗ approximate TW10
10
10
contour plots (bottom) the green squares denote set K.
Table 1: Approximate solutions r̂x∗0 (K, A) are given in the table for several initial conditions
x0 , together with the related sub-optimal action at the initial time.
x0
r̂x∗0 (K, A)
a

[19 19]T
0.8808

[20.5 19]T
0.9454

[19 20.5]T
0.9206

[20.5 20.5]T
0.9557

(ON,ON)

(OFF,ON)

(ON,OFF)

(OFF,OFF)

18]T

21.5]T

[21.5 21.5]T
0.8312
(OFF,OFF)

18]T

x0
[18
∗
r̂x0 (K, A)
0.5204
(ON,ON)
a

[21.5
0.7635
(ON,OFF)

[18
0.8596
(OFF,ON)

6.3 Performance of the Fitted Value Iteration
We are interested in the performance of the FVI algorithm and in analyzing how the computed accuracy deteriorates over the iterations from Nt − 1 to 1. Note that the last iteration
is of little interest, since it does not include the fitting step. The accuracy is computed using
the model-based and sample-based bounds of Section 5. Fig. 3 plots the sample-based esti∗ k
mates of the single step error (19), namely kŴk∗ − T̂Ŵk+1
1,η̃ and of the bias (21), namely
19

Haesaert, Babuska and Abate

µ̂ ∗5( x )

µ̂ ∗1( x )

22

22

20

20

20

18

x2

22

x2

x2

µ̂ ∗9( x )

18
18

20
x1

22

18
18

20
x1

22

18

20
x1

22

Figure 2: The policy µ̂∗ for k = 9, 5, 1 obtained form the same computations as in Fig.
1. The action (ON,ON) is labeled in red, (OFF, ON) is orange, (ON,OFF) is yellow, and
(OFF,OFF) is blue.

5

Estimate

·10−3

4.5
4
3.5
3

0

1

2

3

4

5
k

6

7

8

9

10

∗ k
Figure 3: Sample-based estimation of the single step error kŴk∗ − T̂Ŵk+1
1,η̃ (×) and
a
∗
a
∗
maxa∈A T̂1 Ŵk+1 − T̂2 Ŵk+1 1,η̃ (◦) at each iteration for k = 9 until k = 1. An indepen-

dent set of samples of size Ñ = 4 · 103 and M̃ = 104 has been used.

∗ − T̂a Ŵ ∗
maxa∈A T̂a1 Ŵk+1
2 k+1 1,η̃ . Observe that the values of both (19) and (21) fall in the
interval between 3 × 10−3 and 5 × 10−3 . The bias estimate (21) appears distributed all over
this interval, whereas there is a noticeable trend in the plot of (19), which suggests that
the first iterations can be fitted more easily. In Fig. 4, the accuracy of the FVI algorithm
propagated over the iterations is given, starting from the first iteration kŴ9∗ − TWN∗ t kη
until the last iteration kŴ1∗ − TNt WN∗ t kη . This accuracy is computed using Theorem 12.
The estimates in Fig. 3 are used to compute the estimate of the accuracy kŴk∗ − TWN∗ t kη
and the accuracy ∆ for a given δ∆ . For each iteration step, it can be observed in Fig. 3
that the error caused by estimating the dynamic programming operator T and by fitting
a function is relatively small (< 10−2 ). However, the error grows exponentially over the
whole horizon: as expected, the accuracy of the algorithm depends strongly on B, which
has been computed numerically and amounts to 3.07.

7. Conclusions and Future Work
This article has investigated the performance of a sample-based approximation scheme for
the synthesis of optimal controllers maximizing the probability of the known “reach-avoid”
specification. The approximate computational scheme is based on the Fitted Value Iteration
20

Quantitative Approximations for Reach-Avoid over General Markov Processes

Accuracy ∆
103

100

10−3

0

1

2

3

4

5
k

6

7

8

9

10

Figure 4: The plot shows the accuracy based on the error propagation over the iterations
estimated with an hold out set of size Ñ = 4·103 and M̃ = 104 . The accuracy ∆, marked as
(×), is given for 1−δ∆ = 0.9. In the graph, the estimate of the accuracy, kŴk∗ −TNt −k WN∗ t k,
is given as (◦).
algorithm, which hinges on random sample extractions.
We are interested in the non-trivial extension to continuous control spaces, as well as in
the assessment of the performance of synthesized approximate policies over the concrete
model. Finally, the development of better sampling distributions that minimize the error
propagation can lead to tighter errors, which can be more relevant in practice. To this
end, the optimal sampling distribution should be used to optimize the scaling factors by
resembling more closely the local stochastic kernels.

References
A. Abate, M. Prandini, J. Lygeros, and S. Sastry. Probabilistic reachability and safety
for controlled discrete time stochastic hybrid systems. Automatica, 44(11):2724–2734,
November 2008.
A. Abate, J.P. Katoen, J. Lygeros, and M. Prandini. Approximate model checking of
stochastic hybrid systems. European Journal of Control, 16(6):624–641, 2010.
A. Abate, J.-P. Katoen, and A. Mereacre. Quantitative automata model checking of autonomous stochastic hybrid systems. In Proceedings of the 14th ACM international conference on Hybrid Systems: computation and control, pages 83–92, 2011.
M. Anthony. Uniform Glivenko-Cantelli theorems and concentration of measure in the
mathematical modelling of learning. Research report, Department of Mathematics London School of Economics, May 2002.
M. Anthony and P. L. Bartlett. Neural Network Learning: Theoretical Foundations. cambridge university press, 1999.
C. Baier and J.-P. Katoen. Principles of Model Checking. The MIT Press, 2008.
P.L. Bartlett, O. Bousquet, and S. Mendelson. Local Rademacher Complexities. Annals of
Statistics, 33:1497–1537, 2005.
21

Haesaert, Babuska and Abate

D.P. Bertsekas and S.E. Shreve. Stochastic Optimal Control: The discrete time case. Athena
Scientific, 1996.
H.A.P. Blom and J. Lygeros. Stochastic Hybrid Systems: Theory and Safety Critical Applications. Number 337 in Lecture Notes in Control and Information Sciences. Springer
Verlag, Berlin Heidelberg, 2006.
L. Busoniu, R. Babuska, B.D. Schutter, and D. Ernst. Reinforcement Learning and Dynamic
Programming Using Function Approximators. Automation and Control Engineering. CRC
Press, 2010.
C.G. Cassandras and J. Lygeros. Stochastic Hybrid Systems. Number 24 in Control Engineering. CRC Press, Boca Raton, 2006.
S. Esmaeil Zadeh Soudjani and A. Abate. Adaptive and sequential gridding procedures
for the abstraction and verification of stochastic processes. SIAM Journal on Applied
Dynamical Systems, 12(2):921–956, 2013.
A.M. Farahmand, R. Munos, and C. Szepesvari. Error Propagation for Approximate Policy
and Value Iteration. In Neural Information Processing Systems, pages 568–576, 2010.
A. Fehnker and F. Ivančić. Benchmarks for hybrid systems verification. In R. Alur and
G.J. Pappas, editors, In Hybrid Systems: Computation and Control (HSCC 2004), pages
326–341. Springer, 2004.
E.J. Hartman, J.D. Keeler, and J.M. Kowalski. Layered neural networks with gaussian
hidden units as universal approximations. Neural computation, 2(2):210–215, 1990.
D. Haussler. Decision theoretic generalizations of the PAC model for neural net and other
learning applications. Information and Computation/information and Control, 100:78–
150, 1992.
W. Hoeffding. Probability inequalities for sums of bounded random variables. Journal of
the American Statistical Association, 58(301):13–30, 1963.
M. Karpinski and A. Macintyre. Polynomial bounds for VC dimension of sigmoidal and
general Pfaffian neural networks. Journal of Computer and System Sciences, 54(1):169 –
176, 1997.
M.J. Kearns and R.E. Schapire. Efficient distribution-free learning of probabilistic concepts.
Journal of Computer and System Sciences, 48(3):464 – 497, 1994.
K. Koutsoukos and D. Riley. Computational methods for reachability analysis of stochastic
hybrid systems. In J. Hespanha and A. Tiwari, editors, Hybrid Systems: Computation and
Control, number 3927 in Lecture Notes in Computer Sciences, pages 377–391. SpringerVerlag, Berlin, 2006.
H. J. Kushner and P.G. Dupuis. Numerical Methods for Stochastic Control Problems in
Continuous Time. Springer-Verlag, New York, 2001.
22

Quantitative Approximations for Reach-Avoid over General Markov Processes

S.P. Meyn and R.L. Tweedie. Markov chains and stochastic stability. Springer Verlag, 1993.
R. Munos and C. Szepesvari. Finite time bounds for fitted value iteration. Journal of
Machine Learning Research, 9:815–857, 2008.
Leonid Peshkin and Sayan Mukherjee. Bounds on sample size for policy evaluation in
Markov environments. In Computational Learning Theory, pages 616–629. Springer, 2001.
D. Pollard. Convergence of stochastic processes. Springer Verlag, 1984.
Maria Prandini and Jianghai Hu. Stochastic reachability: Theory and numerical approximation. Stochastic hybrid systems, Automation and Control Engineering Series, 24:107–138,
2006.
S. Summers and J. Lygeros. Verification of discrete time stochastic hybrid systems: A
stochastic reach-avoid decision problem. Automatica, 46:1951–1961, 2010.
I. Tkachev and A. Abate. Characterization and computation of infinite-horizon specifications over markov processes. Theoretical Computer Science, 515:1 – 18, 2014.
V. Vapnik. Statistical Learning Theory. Adaptive and learning systems for signal processing,
communications, and control. John Wiley & Sons, Inc., 1998.

23

Haesaert, Babuska and Abate

Appendix A. Proof of Lemma 6: Bound on Estimation Error
We employ results on the concentration of random variables (Anthony, 2002), which in
general raise conditions on a random variable ensuring its realizations to be concentrated
around its expectation, in the sense that the probability of a deviation from the expectation
is exponentially small (as a function of the deviation). Of interest to this work is a known
bound holding for sums of bounded and independent random variables Hoeffding (1963).
Proposition 13 (Hoeffding’s inequality, Hoeffding (1963)) Suppose that Xi , for i =
1, 2, . . . , N, are independent random variables supported on [0, 1]. Then
o
nP
P
N
i − E N X ≥ N ǫ ≤ 2e−2N ǫ2 ,
P
x
i
i=1
i=1 k

P
P
where E N
Xi is the mean of the random variable i Xi , whereas the empirical mean is
i=1
P
i
i
defined as N
i=1 xk , where xk is a realization of Xi .
Using Proposition 13 the proof of Lemma 6 is provided as follows.

∗
Proof Let us express a probabilistic error bound on the accuracy of the estimate T̂Ŵk+1
i
at each base point xk and given any a ∈ A as
o
n
∗
∗
kp,η̂ > ǫ1 ≤ δ1 ,
− T̂Ŵk+1
P kTŴk+1

where we have used the empirical norm based on η̂. We obtain
o
n
o
n
∗
∗
∗
∗
− T̂Ŵk+1
kpp,η̂ ≤ ǫp1
P kTŴk+1
− T̂Ŵk+1
kp,η̂ ≤ ǫ1 = P kTŴk+1

via definition of the empirical norm in (8)


N
1 P
p
p
∗
i
∗
i
=P
TŴk+1 (xk ) − T̂Ŵk+1 (xk ) ≤ ǫ1
N i=1

Note the mutual independence between the sample sets at different base points xik given as
S
i,a,j 
a∈A xk+1 1≤j≤M
 N n
o Q
n
o
N
p


T
p
∗
i
∗
i
∗
∗ (xi ) ≤ ǫ
TŴk+1 xk − T̂Ŵk+1 (xk ) ≤ ǫ1
≥P
=
P TŴk+1
xik − T̂Ŵk+1
1 .
k
i=1

i=1

Let us now express the argument of the probability operator as follows

∗
∗
TŴk+1
xik − T̂Ŵk+1
(xik )
h
i
∗
= max Exk+1 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
(xk+1 )
a∈A

− max
a∈A

M
1 P
i,a,j
i,a,j
∗
[1K (xi,a,j
k+1 ) + 1A\K (xk+1 )Ŵk+1 (xk+1 )]
M j=1

W.r.t Exk+1 defined over random variable xk+1 ∼ Tx · | xik , a
h
i
∗
≤ max Exk+1 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
(xk+1 )
a∈A

−

M
1 P
i,a,j
i,a,j
∗
[1K (xi,a,j
k+1 ) + 1A\K (xk+1 )Ŵk+1 (xk+1 )] .
M j=1

24



Quantitative Approximations for Reach-Avoid over General Markov Processes

Therefore the probability of the last event above can be lower bounded by the probability
associated to several independent events over the finite action space, as follows:
n
o

∗
∗
P TŴk+1
xik − T̂Ŵk+1
(xik ) ≤ ǫ1
i
h
Q n
∗ (x
P Exk+1 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
≥
k+1 )
a∈A

−

i
o
M h
1 P
i,a,j
∗ (xi,a,j )
.
1K (xi,a,j
)
+
1
(x
)
Ŵ
≤
ǫ
1
A\K
k+1 k+1
k+1
k+1
M j=1

∗
∈ W, define random
For a given base point xik ∈ X , action a ∈ A, and function Ŵk+1
i,a,j
i,a,j
i,a,j
∗
variables Zj via their realizations 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1 (xk+1 ), with j = 1, . . . , M .

i
Since each xi,a,j
k+1 is independently drawn from Tx · | xk , a , the random variables Zj are
independent, identically distributed, and take values within the closed interval [0, 1]. By
application of Hoeffding’s inequality (as in Proposition 13), the concentration of the M
samples around the expected value of Zj can be expressed as
n
h
i
∗
P Exk+1 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
(xk+1 )
i o
M h
1 P
i,a,j
∗ (xi,a,j ) ǫ
−2M (ǫ1 )2 .
−
1K (xi,a,j
)
+
1
(x
)
Ŵ
1 ≤ 2e
A\K k+1
k+1 k+1
k+1
M j=1
o
n
2
∗ k
∗
≤
ǫ
− T̂Ŵk+1
Therefore as long as 0 ≤ 2e−2M (ǫ1 ) ≤ 1, it follows that P kTŴk+1
1 ≥
p,η̂

2 N |A|
.
1 − 2e−2M (ǫ1 )

Remark 14 As long as we only know that the random variables Zj are bounded, the use
of Hoeffding’s inequality is sufficient. If we further have information on the variance of
Zj , one can leverage the inequalities of Chebyshev and of Bienaym-Chebyshev (Hoeffding,
1963), or alternatively Bernstein’s inequality (Peshkin and Mukherjee, 2001): the former
bounds are only function of the variance, whereas the latter inequality depends not only on
variance of Zj but also on its bounded domain. Upper bounds on either the variance of Zj
or on its range can be derived exploiting prior knowledge on properties of the function space
W and of the distribution Tx (· | x, a).

Appendix B. Proof of Lemma 7
We derive a general, analytical bound on the error of a single backward recursion using
notions from statistical learning theory (Haussler, 1992; Pollard, 1984). The error bound
∗ , the optimal fit can be anywhere in the function
takes into account that, for any TŴk+1
class. Furthermore the bound will be distribution-free, namely holding for any Markov
process (with dynamics characterized by Tx ) and any sample distribution η over the set
A \ K.
25

Haesaert, Babuska and Abate

We exclusively consider function classes W ⊂ B(X ; 1) endowed with a finite pseudodimension: this includes all finitely-parameterized function classes (Munos and Szepesvari,
2008). The notion of pseudo dimension (Pollard, 1984; Anthony, 2002; Haussler, 1992)
expresses the capability of a function class W to the fit a set of samples.
Proof In order to prove Lemma 7, we show that the inequality in (10) holds for any
∗
∈ W at any time instant k = 0, . . . , Nt − 1. For the sake of notation in the following
Ŵk+1
∗
we substitute Ŵk+1
by W , and instead of considering the set of base points (xik )1≤i≤N drawn
at the time instant k we simply introduce ~x = (x1 , . . . , xN ) as a sequence of N independent
realizations drawn from a distribution over A \ K with density η.
For any given function W ∈ W , induce a new function class lW = {|w − TW |p : w ∈ W}
with elements lw ∈ lW : lw = |w − TW |p . The inequality in (10) can be rewritten over the
function class lW as follows
o
n
o
n
N
1 P
lw (xi ) ≥ ǫp2 ,
P sup kw − TW kpp,η − kw − TW kpp,η̂ ≥ ǫp2 = P sup Eη [lw ] −
N i=1
w∈W
lw ∈lW
where Eη denotes the expected value with respect to η. This allows us to use a result in
(Pollard, 1984), which provides an upper-bound on the probability of the above event as a
function of the covering number of the metric space ((lW )~x , k · k1 ).
Proposition 15 (Pollard (1984)) Let F be a permissable set of functions on X with
0 ≤ f (x) ≤ K for all f ∈ F and x ∈ X . Let ~x = (x1 , . . . , xN ) be a sequence of N samples
drawn independently from X according to any distribution on X . Then for all ǫ > 0
n
P ∀f ∈ F : Ef −

1
N

Nǫ2

−
2
x , k · k1 )e 128K ,
~
x f (xi ) ≥ ǫ ≤ 4EN (ǫ/16, F|~

P

(26)

where the quantity N will be introduced shortly and where the definition of a permissible set
of functions (Pollard, 1984) includes all finitely parameterized functions.
Let us introduce the concept of covering number of a metric space (Haussler, 1992).
Given a (pseudo-)metric space (A, ρ) and a subset S of A, we say that the set T ⊆ A is an
ǫ-cover for S (where ǫ > 0) if, for every s ∈ S there is a t ∈ T such that ρ(s, t) < ǫ. For a
given ǫ > 0 we denote the covering number N (ǫ, S, ρ) (Haussler, 1992) as the cardinality of
the smallest ǫ-cover of S.
For a given set of samples xi with i = 1, . . . , N , the evaluation of a function lw ∈
lW over each of these samples is given as the N dimensional vector in [0, 1]N : (lw )|~x =
(lw (x1 ), lw (x2 ), . . . , lw (xN )). The induced set of vectors is
(lW )|~x = {(lw )|~x = (lw (x1 ), lw (x2 ), . . . , lw (xN )), lw ∈ lW } ⊆ [0, 1]N .
The minimal ǫ-cover of ((lW )|~x , k · k1 ) is denoted as N (ǫ, (lW )|~x , k · k1 ).
The deviation of the expected value from the empirical mean can be bounded using
Pollard’s proposition (Pollard, 1984)
n
P suplw ∈lW Ex [lw (x)] −

1
N

o
h
i N(ǫ2 )2p
p
p
− 128
l
(x
)
≥
ǫ
≤
4E
N
(ǫ
/16,
(l
)
,
k
·
k
)
.
W |~
1 e
2
2
i=1 w i
x

PN

26

Quantitative Approximations for Reach-Avoid over General Markov Processes

The expected value of N (ǫp2 /16, (lW )|~x , k · k1 ) is computed over the samples xi of ~x, drawn
independently from a probability distribution with density η. Since there is a trivial isometry
(Haussler, 1992) between (lW |~x , k · k1 ) and (lW , k · k1,η̂ ), both spaces have equal covering
numbers
N (ǫp2 /16, lW |~x , k · k1 ) = N (ǫp2 /16, lW , k · k1,η̂ ).
In practice a value for E [N (ǫp2 /16, lW , k · k1,η̂ )] can be obtained by upper bounding
N (ǫp2 /16, lW , k · k1,η̂ ) independently of the sample distribution. For this we introduce the
pseudo dimension of a function class, formally defined as follows Pollard (1984); Anthony
(2002); Haussler (1992). Suppose F is a class of functions, f ∈ F, f : X → [0, 1]. Then
S ⊆ X is shattered by F if there are numbers rx ∈ [0, 1] for x ∈ S such that for every T ⊆ S
there is some fT ∈ F with the property that fT ≥ rx if x ∈ T and fT < rx if x ∈ S \ T . We
say that F has a finite pseudo dimension dimp (F) = d if d is the maximum cardinality of
a shattered set.
For any distribution P ∈ M (X ), the packing number (Haussler, 1992) and therefore also
tho covering number of the metric space (lW , k · k1,P ) can be upper bounded as a function
of the pseudo-dimension and the base of the natural logarithm e: for any ǫ > 0,
N (ǫ, lW , k · k1,P ) ≤ e(d + 1)


2e d
,
ǫ

with dimp (lW ) = d.

We have proved that a sufficient upper bound is given as
n
o
 d N(ǫ2 )2p
e− 128 .
P supw∈W kw − TW kpp,η − kw − TW kpp,η̂ ≥ ǫp2 ≤ 4e(d + 1) 32e
ǫp
2

The proof can be concluded by showing that the pseudo dimension d of the induced class lW
is the same as the pseudo dimension of W. Let {w − TW : w ∈ W} be a new function class
induced from W. The invariance propertiesof the pseudo dimension dimp (W) shown in
(Haussler, 1992) allow to conclude that dimp ( w − TW w ∈W ) = dimp (W). The induced
function class lW can then be defined as follows: lW = |k|p k ∈ {w − TW : w ∈ W} .
Since it was shown in (Kearns and Schapire, 1994) that the pseudo dimension is invariant
over function composition (|·|p ), we conclude that the pseudo dimension is dimp (lW ) =
dimp ({w − TW : w ∈ W}) = dimp (W) = d.
Remark 16 (Computing the pseudo-dimension) When the function class W is a vector space of real-valued functions, the pseudo dimension is equal to the dimensionality of the
function class (Anthony and Bartlett, 1999, Theorem 11.4). (Anthony and Bartlett, 1999)
elaborates the details of the computation of pseudo dimensions of parameterized function
classes, especially for function classes defined over neural networks.
Since it is possible to bound the pseudo dimension of lW (as introduced in the proof) by
the pseudo dimension of W, this capacity concept has been used to bound the error caused
by using an empirical estimate of the weighted p-norm. Notice that for non-parametric
function classes, concepts such as covering number or Rademacher average of the function
class lW can be used instead (Bartlett et al., 2005).
27

Haesaert, Babuska and Abate

Let us shortly discuss how the derived bounds can be tightened. A first option is to
circumvent the notion of pseudo dimension and work with the covering numbers in Pollard
inequality (Proposition 15), however the increase in assumptions on the function class and
in overall computations make the gain in accuracy undeserving. A second option is to
explore alternatives over Pollard inequality in (15) with better constants (Bartlett et al.,
2005). An alternative concentration inequality based on Bernstein’s inequality is used in
(Peshkin and Mukherjee, 2001). Hoeffding inequality gives a concentration inequality on
the sum of bounded random variables, whereas Bernstein inequality gives a tighter bound
based on knowledge of both the boundness and the variance of the random variables. Even
with improved constants or alternative inequalities, the error bounds can still result to be
conservative for reasonable sample complexities.

Appendix C. Proof of Theorem 8
The proof of Theorem 8 is adapted from the proof of the single-step error bound for Fitted
Value Iteration with multiple sample batches in (Munos and Szepesvari, 2008).
∗
Proof Let us introduce a simplified notation for Ŵk+1
by replacing it with a general
′
′
function W ∈ W that minimizes the empirical norm as W = arg minw∈W kw− T̂W kp,η̂ . Let
us further define a space Ω for the batch of samples drawn atany ofthe iterations, such that



S
S
i,a,j
i
at any instant k the realized sample batch ω := i∈{1,...,N } xk ∪
a∈A xk+1
1≤j≤M

is an element of the sample space, ω ∈ Ω.
For any given ǫ′ > 0, consider a function w∗ ∈ W such that kw∗ −TW kp,η ≤ inf w∈W kw−
TW kp,η + ǫ′ (this in particular holds since W has been assumed to be close and bounded).
The error bound in (13) holds for a sample realization ω if the following sequence of
inequalities holds simultaneously:
kW ′ − TW kp,η ≤ kW ′ − TW kp,η̂ + ǫ2

≤ kW ′ − T̂W kp,η̂ + ǫ1 + ǫ2
≤ kw∗ − T̂W kp,η̂ + ǫ1 + ǫ2

(27a)
(27b)
(27c)

∗

(27d)

∗

(27e)

≤ kw − TW kp,η̂ + 2ǫ1 + ǫ2

≤ kw − TW kp,η + 2ǫ1 + 2ǫ2 .

As long as the previous sequence of inequalities is true, the following one also holds:
kW ′ − TW kp,η ≤ inf kw − TW kp,η + 2ǫ1 + 2ǫ2 + ǫ′ .
w∈W

We claim that the sequence of inequalities holds with a probability at least 1−(δ1 +δ2 ). Since
there exists a function w∗ for any ǫ′ > 0 it follows with a probability at least 1 − (δ1 + δ2 )
that
kW ′ − TW kp,η ≤ dp,η (TW, W) + 2ǫ1 + 2ǫ2 .
By the union bound argument (Anthony, 2002), the probability of the union of events can
be bounded by the sum of the probabilities of the single events. Using this argument it
is possible to define a lower bound on the probability associated with the simultaneous
28

Quantitative Approximations for Reach-Avoid over General Markov Processes

occurrence of the five inequalities in (27). We first show that the third inequality is always
true. Then we give the probability associated to the first inequality (27a) and the fifth
(27e) (this is based on (12)). Afterwards we provide an upper bound on the probability
associated to the second and fourth inequalities (27b),(27d), based on the bound given in
(11).
The third inequality (27c) is true for the whole sample space Ω due to the choice of
W ′ . For all functions w in W it follows that kW ′ − T̂W kp,η̂ ≤ kw − T̂W kp,η̂ holds, because
W ′ = arg minw∈W kw − T̂W kp,η̂ .
The first and last inequalities (27a),(27e) bound the deviation between the empirical
loss and the expected loss. This can be bounded with the worst case error. Firstly we
observe that the inequality
p

kw − TW kp,η̂ − kw − TW kp,η

≤ kw − TW kpp,η − kw − TW kpp,η̂

is always true. In the case that kw − TW kp,η̂ ≤ kw − TW kp,η then
kw − TW kp,η̂ − kw − TW kp,η = kw − TW kp,η̂ − kw − TW kp,η
kw − TW kp,η̂ = (kw − TW kp,η̂ − kw − TW kp,η ) + kw − TW kp,η

kw − TW kpp,η̂ = ((kw − TW kp,η̂ − kw − TW kp,η ) + kw − TW kp,η )p

kw − TW kpp,η̂ ≥ (kw − TW kp,η̂ − kw − TW kp,η )p + kw − TW kpp,η

kw − TW kpp,η̂ − kw − TW kpp,η ≥ (kw − TW kp,η̂ − kw − TW kp,η )p
kw − TW kpp,η̂ − kw − TW kpp,η ≥ kw − TW kp,η̂ − kw − TW kp,η

p

On the other hand, for the case when kw − TW kp,η̂ > kw − TW kp,η a similar argument can
be used. We can then observe that
kw∗ − TW kp,η̂ − kw∗ − TW kp,η

p

≤ sup kw − TW kpp,η − kw − TW kpp,η̂ ,
w∈W

and that
kW ′ − TW kp,η − kW ′ − TW kp,η̂

p

≤ sup kw − TW kpp,η − kw − TW kpp,η̂ .
w∈W

Given two functions w∗ and W ′ define events A1 and A2
p

A1 : ǫp2 < kw∗ − TW kp,η̂ − kw∗ − TW kp,η ,

Observe that the event sets A1 and A2 are subsets of the more general event B defined as
B : ǫp2 < sup kw − TW kpp,η − kw − TW kpp,η̂ .
w∈W

Thus it follows that for any ǫ2 > 0: P {A1 ∪ A2 } ≤ P {B} and, based on (12), we have
nn
o n
oo
P
kW ′ − TW kp,η − kW ′ − TW kp,η̂ > ǫ2 ∪ kw∗ − TW kp,η̂ − kw∗ − TW kp,η > ǫ2


p
p
p
≤ P sup kw − TW kp,η − kw − TW kp,η̂ > ǫ2 ≤ δ2 .
w∈W

29

p

A2 : ǫp2 < kW ′ − TW kp,η − kW ′ − TW kp,η̂ .

Haesaert, Babuska and Abate

Thus the probability that the inequalities (27a) and (27e) do not hold is less then δ2 .
The second and fourth inequalities (27b),(27d) depend the accuracy of the estimation of
the backward recursion at each base point xik . Employing the inequality |kw − gkp,η̂ − kw − hkp,η̂ | ≤
kg − hkp,η̂ , we can see that
kW ′ − TW kp,η̂ − kW ′ − T̂W kp,η̂ ≤ kTW − T̂W kp,η̂ ,
and
kw∗ − T̂W kp,η̂ − kw∗ − TW kp,η̂ ≤ kTW − T̂W kp,η̂ .
For every sample set ω the inequalities (27b),(27d) apply if kTW − T̂W kp,η̂ ≤ ǫ1 . Thus
P

nn

oo
o n
kW ′ − TW kp,η̂ − kW ′ − T̂W kp,η̂ > ǫ1 ∪ kw∗ − T̂W kp,η̂ − kw∗ − TW kp,η̂ > ǫ1
o
n
≤ P kTW − T̂W kp,η̂ > ǫ1 ≤ δ1 .

(28)

The probability that at least one of the inequalities in (27) does not hold can be expressed
using the union bound as δ1 + δ2 . Thus the sequence of inequalities holds with at least a
probability of 1 − δ1 − δ2 .

Appendix D. Proof of Lemma 9
Proof Let us set up the following chain of inequalities:
kTNt −k ŴN∗ t − Ŵk∗ kp,η

∗
]
= [ Add and subtract function TŴk+1


∗
∗
− Ŵk∗ kp,η
+ Ŵk+1
= kT TNt −k−1 ŴN∗ t − Ŵk+1

= [ Definition of T in (2), where we have considered a single xk ∼ η ]

i
h

∗
∗
(xk+1 ) xk+1 ∼ Tx (· | xk , a)
+ Ŵk+1
= max E 1K (xk+1 ) + 1A\K (xk+1 ) TNt −k−1 ŴN∗ t − Ŵk+1
a∈A

− Ŵk∗

p,η

= [ max E[ξ1 + ξ2 ] ≤ max E|ξ1 | + max E|ξ2 | ]


∗
≤ max E 1A\K (xk+1 ) TNt −k−1 ŴN∗ t − Ŵk+1
(xk+1 )|xk+1 ∼ Tx (· | xk , a)
a∈A

∗
(xk+1 )|xk+1 ∼ Tx (· | xk , a) − Ŵk∗
+ max E 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
a∈A

= [ Triangular inequality]


∗
≤ max E 1A\K (xk+1 ) TNt −k−1 ŴN∗ t − Ŵk+1
(xk+1 )|xk+1 ∼ Tx (· | xk , a)
a∈A

30

p,η

p,η

Quantitative Approximations for Reach-Avoid over General Markov Processes

∗
+ max E 1K (xk+1 ) + 1A\K (xk+1 )Ŵk+1
(xk+1 )|xk+1 ∼ Tx (· | xk , a) − Ŵk∗
a∈A

= [ Definition of T in (2) ]

i
h

∗
= max E 1A\K (xk+1 ) TNt −k−1 ŴN∗ t − Ŵk+1
(xk+1 ) |xk+1 ∼ Tx (· | xk , a)
a∈A

∗
− Ŵk∗
+ TŴk+1

p,η

p,η

= [ Introduce density function tx (xk+1 | xk , a) for kernel Tx ]
Z


∗
= max
(xk+1 ) tx (xk+1 |xk , a)dxk+1
1A\K (xk+1 ) TNt −k−1 ŴN∗ t − Ŵk+1
a∈A

p,η

X

p,η

1

∗
Let us now show that the first term is bounded by B p TNt −k−1 ŴN∗ t − Ŵk+1

max
a∈A

Z

X



∗
tx (xk+1 |xk , a)dxk+1
1A\K (xk+1 ) TNt −k−1 ŴN∗ t − Ŵk+1

X

= [ Express the η-weighted p-norm over A \ K ]

p
R
R
∗
tx (xk+1 |xk , a)dxk+1
= A\K maxa∈A A\K TNt −k−1 ŴN∗ t − Ŵk+1
=



R

A\K

maxa∈A

R

Nt −k−1 Ŵ ∗
Nt
A\K T

−

∗
Ŵk+1

p

p,η

:

p,η

= [ monotonicity of Lp -norms with respect to a probability measure ]
Z
1
p

p
Nt −k−1
∗
∗
≤ max
1A\K (xk+1 ) T
tx (xk+1 |xk , a)dxk+1
ŴNt − Ŵk+1
a∈A

∗
TŴk+1
− Ŵk∗

tx (xk+1 |xk , a)dxk+1

1
p

p

1
p

p

p,η

η(xk )dxk
η(xk )dxk

1
p
p
Nt −k−1 Ŵ ∗ − Ŵ ∗
T
t
(x
|x
,
a)dx
η(x
)dx
x k+1 k
k+1
k
k
Nt
k+1
A\K

R

1
p
R
p
∗
≤ A\K A\K maxa∈A TNt −k−1 ŴN∗ t − Ŵk+1
tx (xk+1 |xk , a) dxk+1 η(xk )dxk
=

R

A\K maxa∈A

R

η(x

 p1
 p1

)

= [ Introduce dummy term η(xk+1
, which is defined over xk+1 ∈ A \ K ]
k+1 )
1
R


p
R
p
η(xk )
∗
dx
η(x
)dx
= A\K A\K maxa∈A TNt −k−1 ŴN∗ t − Ŵk+1
tx (xk+1 |xk , a) η(x
k
k+1
k+1
k+1 )
p

∗
is only a function of xk+1 ]
= [ Recall that TNt −k−1 ŴN∗ t − Ŵk+1

1

R
pR
p
tx (xk+1 |xk ,a)η(xk )
∗
.
dx
η(x
)dx
max
= A\K TNt −k−1 ŴN∗ t − Ŵk+1
a∈A
k
k+1
k+1
A\K
η(xk+1 )



R
t (x
|xk ,a)η(xk )
dxk over the domain
Introduce now the upper bound on A\K maxa∈A x k+1
η(xk+1 )
R
t (x
|xk ,a)η(xk )
dxk , obtaining
A \ K as B = supxk+1 ∈A\K A\K maxa∈A x k+1
η(xk+1 )
≤

R

A\K

p

∗
TNt −k−1 ŴN∗ t − Ŵk+1
Bη(xk+1 )dxk+1

31

1

p

1

∗ k .
= B p kTNt −k−1 ŴN∗ t − Ŵk+1
p,η

p,η

.

Haesaert, Babuska and Abate

We have finally shown that TNt −k ŴN∗ t − Ŵk∗

p,η

∗
− Ŵk∗
≤ TŴk+1

1

p,η

∗
+B p TNt −k−1 ŴN∗ t − Ŵk+1

Appendix E. Proof of Theorem 10



Proof If we estimate the quantity rx∗0 (K, A) = TNt WN∗ t (x0 ) = TNt ŴN∗ t (x0 ) by




T̂Ŵ1∗ (x0 ), then we have that r̂x∗0 (K, A) = 1K (x0 ) + 1A\K (x0 ) T̂Ŵ1∗ (x0 ). The absolute
deviation of the approximated r̂x∗0 (K, A) from the exact rx∗0 (K, A) is given as




r̂x∗0 (K, A) − rx∗0 (K, A) = TNt ŴN∗ t (x0 ) − T̂Ŵ1∗ (x0 ) .

The objective is to 
present this
function of the errors introduced by the
 erroras a 
N
∗
∗
t
approximate mappings T ŴNt (x0 ) − T̂Ŵ1 (x0 ) , as well as of the quantities kŴ1∗ −
TŴ2∗ kp,η , kŴ2∗ − TŴ3∗ kp,η , . . . , kŴN∗ t −1 − TŴN∗ t kp,η .




To this end, we first express a bound on TNt ŴN∗ t (x0 ) − T̂Ŵ1∗ (x0 ) as a function


 
. Then Lemma 9 is used to
of | T̂ TŴ1∗ (x0 ) − Ŵ1∗ (x0 )| and of TNt −1 ŴN∗ t − Ŵ1∗
p,η

express kT Nt −1 ŴN∗ t − Ŵ1∗ kp,η as a function of the errors introduced by the approximate
mappings. Similar to the first chain of inequality in the proof of Lemma 9 applied at step
k = 0 and point x0 , we obtain that



TNt ŴN∗ t (x0 ) − T̂Ŵ1∗ (x0 )
Z




≤ max
TNt −1 ŴN∗ t (x1 ) − Ŵ1∗ (x1 ) tx (x1 | x0 , a) dx1 + TŴ1∗ (x0 ) − T̂Ŵ1∗ (x0 ) .


a∈A

A\K

Let us now introduce a measure for the maximum concentration of the density function
tx (x1 | x0 , a) over x1 ∈ A \ K, for any a ∈ A, defined relative to the density of the distri1 |x0 ,a)
. Since B0 η(x1 ) ≥ tx (x1 | x0 , a), it
bution η in (16), as B0 = supx1 ∈A\K maxa∈A tx (xη(x
1)
follows that
max
a∈A

Z

A\K

TNt −1 ŴN∗ t (x1 ) − Ŵ1∗ (x1 ) tx (x1 | x0 , a) dx1 ≤ B0

Z

A\K

TNt −1 ŴN∗ t (x1 ) − Ŵ1∗ (x1 ) η(x1 )dx1 .

The last expression corresponds to a 1-norm with respect to a probability measure η over
A \ K. Exploiting the monotonicity of the p-norm with respect to a probability measure, a
more general expression for the approximation error is obtained as








.
TNt ŴN∗ t (x0 ) − T̂Ŵ1∗ (x0 ) ≤ T̂Ŵ1∗ (x0 ) − TŴ1∗ (x0 ) + B0 TNt −1 ŴN∗ t − Ŵ1∗
p,η

The second term can be expressed as a function of the weighted p-norm of the approximations by applying Lemma 9. This leads to the expression for an upper bound on the
approximation error as




PNt −1 k−1
∗
.
B p Ŵk∗ − TŴk+1
r̂x∗0 (K, A) − rx∗0 (K, A) ≤ T̂Ŵ1∗ (x0 ) − TŴ1∗ (x0 ) + B0 k=1
p,η

32

p,η

.

Quantitative Approximations for Reach-Avoid over General Markov Processes

From the above expression, a sufficient condition the accuracy in (3) to hold is
N
t −1



n 
X
k−1
∗
B p Ŵk∗ − TŴk+1
P
T̂Ŵ1∗ (x0 ) − TŴ1∗ (x0 ) + B0
k=1

p,η

o
> ∆ ≤ δ∆ .

Appendix F. Sample Complexities
Given ǫ0,1,2 and α, select δ0,1,2 > 0 such that 1 − α = δ0 + (Nt − 1)δ1 + (Nt − 1)δ2 , and let
us pick values for N ,M , M0 such that
 32e d Nǫ2p
2
−2M0 (ǫ0 )2
−2M (ǫ1 )2
δ0 ≤ 2|A|e
, δ1 ≤ 2|A|N e
, δ2 ≤ 4e(d + 1) p
e− 128 .
ǫ2

Note that the first two inequalities are approximated with first order approximation for
2
2
2
which we know that 1 − (1 − 2e−2M0 (ǫ0 ) )|A| ≤ 2|A|e−2M0 (ǫ0 ) and 1 − (1 − 2e−2M (ǫ1 ) )|A|N ≤
2
2|A|N e−2M (ǫ1 ) . The obtained integer values for N ,M , M0 are given as

 2p   m
 2p  
l
 2p

1
1
1
1

+
128
+
128dp
ln
ln δ12 ,
N
=
128
(ln(4e(d
+
1))
+
d
ln(32e))

ǫ2
ǫ2
ǫ2
ǫ2


m
l  2 
1
1
.
M =
ln(2|A|) + ln( δ11 ) + ln(N ) ,
2 ǫ1


l  2 
m


1
1
 M0 =
ln(2|A|) + ln( δ10 ) ,
2 ǫ0

The use of the obtained M, M0 , N in (17) leads to a confidence of at least α.

Appendix G. Proof of Theorem 12
Proof The proof of Theorem 12 is built observing that (a.) the single step error kŴk∗ −
∗ k
TŴk+1
1,η is bounded by the sum of the expectations of (19) and (21); that (b.) the
propagation of the single step errors gives a bound on the overall approximation error, see
Theorem 10 – hence the expected value of the estimates, propagated over the time horizon,
also gives a bound on the approximation error; and that (c.) the one-sided application
of the Hoeffding’s inequality provides a probabilistic upper bound on the deviation of the
estimate from its mean, and therefore also bounds the approximation error probabilistically.
Part (a.)
i
h
∗
∗
with Ex [f (x)] the mean of f (x) for x ∼ η.
(x)
kŴk∗ − TŴk+1
k1,η = Ex Ŵk∗ (x) − TŴk+1

a,1 a,2
a,M̃
Define a set of i.i.d. random variables
h ~y1 = [y1 , y1 , . . . ,iy1 ] drawn from the distribution
∗ (x)|x as an auxiliary variable with T̂a
y1a,j ∼ Tx (· | x, a). Introduce E~y1 maxa∈A T̂a1 Ŵk+1
1
the estimated operator as defined in (20) and computed over the ~y1 .
= Ex








∗
∗
∗
a
Ŵk∗ (x) − E~y1 max T̂a
y1 max T̂1 Ŵk+1 (x)|x − TŴk+1 (x)
1 Ŵk+1 (x)|x + E~
a∈A

a∈A

33

Haesaert, Babuska and Abate

≤ Ex









∗
a
∗
∗
Ŵk∗ (x) − E~y1 max T̂a
Ŵ
(x)|x
+
E
E
T̂
Ŵ
(x)|x
−
T
Ŵ
(x)
max
x
y
~
1 k+1
1 k+1
k+1
1

≤ Ex,y~1
|

a∈A



a∈A



∗
Ŵk∗ (x) − max T̂a
1 Ŵk+1 (x) + Ex
a∈A
{z
} |
[ Single step error ]

Observe that the single step error, Ex,y~1

h







∗
∗
E~y1 max T̂a
1 Ŵk+1 (x)|x − TŴk+1 (x) .
a∈A
{z
}
[ Bias term ]

∗ (x)
Ŵk∗ (x) − maxa∈A T̂a1 Ŵk+1

i

is equal to kŴk∗ −

∗
∗
∗ k
T̂Ŵk+1
1,η and EkŴk − T̂Ŵk+1 k1,η̃ . The bias term gives the bias introduced by using an
estimate of the operator and it can be rewritten as the expected value of (21). Note that
maxa∈A Ey [Vk+1 (y)|x, a] is a function of x, a and |E~y1 [f (~y )] | ≤ E~y1 [|f (~y )|], thus it follows
that
[ Bias term ] ≤ Ex E~y1



∗
∗
max T̂a
1 Ŵk+1 (x) − TŴk+1 (x)
a∈A



i
h
∗
a
∗
x ≤ Ex E~y1 max T̂a
~2 T̂2 Ŵk+1 (x)|x
1 Ŵk+1 (x) − Ey


x .

a∈A

The second inequality follows from introducing a secondary set of random variables ~y2 =
[y2a,1 , y2a,2 , . . . , y2a,1M̃ ] for which the elements are i.i.d. as y2a,1h ∼ Tx (· | x, a) and
i which are
∗
a
∗
independent of ~y1 . Substituting TŴk+1 (x) with maxa∈A E~y2 T̂2 Ŵk+1 (x)|x, a we have

i 
h
∗
a
∗
[ Bias term ] ≤ Ex E~y1 max E~y2 T̂a
y1 |x
1 Ŵk+1 (x) − T̂2 Ŵk+1 (x) |x, a, ~
a∈A


∗
a
∗
= Ex,~y1 ,~y2 max T̂a
1 Ŵk+1 (x) − T̂2 Ŵk+1 (x) .
a∈A

∗ −
The last equality is equal to the expected value of the estimated bias E maxa∈A T̂a1 Ŵk+1
∗
T̂a2 Ŵk+1
. This proves statement (a.).
1,η̃

if

Part (b.) & (c.) Based on Theorem 10 r̂x∗0 (K, A) has accuracy ∆ with probability δ∆



NX
t −1
n
o
∗
a
∗
a
∗
B k−1 E Ŵk∗ (x) − max T̂a
P |Ŵ0∗ (x0 ) − TŴ1∗ (x0 )| + B0
≥ ∆ < δ∆ .
1 Ŵk+1 (x) + max T̂1 Ŵk+1 (x) − T̂2 Ŵk+1 (x)
k=1

a∈A

a∈A

Which holds under a union bounding argument if
n
o
P |Ŵ0∗ (x0 ) − TŴ1∗ (x0 )| ≥ ǫ0 < δ0
(29)


i
h
P
Nt −1 k−1
a
∗
∗

Ŵ ∗ (x) + maxa∈A T̂a
≥ B0 ǫ 
B
E Ŵk∗ (x) − maxa∈A T̂a

 B0 k=1
2
1 Ŵk+1 (x) − T̂2 Ŵk+1 (x)
−2 Ñ ǫ2

 1 k+1
L
<e
P
P
Nt −1 k−1
∗
∗
a
∗
∗
a


+ maxa∈A T̂1 Ŵk+1 − T̂2 Ŵk+1
Ŵk − T̂Ŵk+1
+B0 k=1 B


1,η̃

1,η̃

(30)

and ∆ and δ∆ are given as (22a) and (22b).
The probabilistic bound (29) follows from Lemma 8 for the estimation error of an empirical
2
norm with accuracy ǫ0 , δ0 obtained for p = 1, M = M0 and N = 1 as long as 0 < 2e−2M0 ǫ0 <
1. The probabilistic bound (30) follows from a one-sided Hoeffding’s inequality (Hoeffding,
1963) with random variable
PNt −1
k=1

B k−1



a
∗
a
∗
∗ (x) + max
Ŵk∗ (x) − maxa∈A T̂a1 Ŵk+1
a∈A T̂1 Ŵk+1 (x) − T̂2 Ŵk+1 (x)

34



,

Quantitative Approximations for Reach-Avoid over General Markov Processes

obtained from the combination of random variable
∼ η and conditional random variables
PNxt −1
~y1 and ~y2 and taking values in the range [0, 2 k=1
B k−1 ]. Note that its estimated of
interest over Ñ samples can be rewritten in the form of (22a),

PNt −1 k−1 
a Ŵ ∗
a Ŵ ∗
∗
−
T̂
T̂
+
Ŵk∗ − T̂Ŵk+1
max
.
a∈A
2
1
k=1 B
k+1
k+1
1,η̃
1,η̃
This concludes the proof of Theorem 12.

Appendix H. Scaling factor for case study
Compute B as in (14) using the given density distribution of the transitions (25), as



Z
exp − 21 (y − µ)T Σ−1 (y − µ) η(x)
1
p
B = sup
max
dx
η(y)
|Σ|(2π)2 a∈A
y∈A\K A\K
= [ µ is a function of a and x, and η(·) is constant over A \ K ]


Z
1
1
T −1
p
max exp − (y − µ) Σ (y − µ) dx
= sup
2
|Σ|(2π)2 a∈A
y∈A\K A\K

= [Suppose A is invertible, and define µ̄(y, a) = A−1 y − A−1 Ba − A−1 C, Σ̄ = A−1 ΣA−T ]


Z
1
1
1
T −1
p
max exp − (x − µ̄(y, a)) Σ̄ (x − µ̄(y, a)) dx
= sup
2
|Σ̄|(2π)2 a∈A
y∈A\K |A| A\K



Z
X
1
1
1
T −1
p
≤ sup
dx
exp − (x − µ̄(y, a)) Σ̄ (x − µ̄(y, a))
2
|Σ̄|(2π)2 a∈A
y∈A\K |A| A\K
 !

XZ
1
1
1
T −1
p
exp − (x − µ̄(y, a)) Σ̄ (x − µ̄(y, a)) dx .
= sup
2
|Σ̄|(2π)2
y∈A\K |A| a∈A A\K
The integral is rewritten as one over a scaled 2-dimension multivariate Gaussian density
distribution with mean µ̄ and covariance Σ̄. With this result, it can be deduced that B is
1
|A| as
smaller than |A|

 
1  P R
√ 1 2 exp − 12 (x − µ̄(y, a))T Σ̄−1 (x − µ̄(y, a)) dx
|Σ̄|(2π)
y∈A\K |A| a∈A X
1 P
1
1 = |A|
|A|.
≤ sup
|A|
a∈A
y∈A\K

B ≤ sup

35

(31)

