text: Optimistic Optimization for Statistical Model Checking with Regret Bounds | DeepAI
  We explore application of multi-armed bandit algorithms to statistical model checking (SMC) of Markov chains initialized to a set of states. We observe that model checking problems requiring maximization of probabilities of sets of execution over all choices of the initial states, can be formulated as a multi-armed bandit problem, for appropriate costs and rewards. Therefore, the problem can be solved using multi-fidelity hierarchical optimistic optimization (MFHOO). Bandit algorithms, and MFHOO in particular, give (regret) bounds on the sample efficiency which rely on the smoothness and the near-optimality dimension of the objective function, and are a new addition to the existing types of bounds in the SMC literature. We present a new SMC tool—HooVer—built on these principles and our experiments suggest that: Compared with exact probabilistic model checking tools like Storm, HooVer scales better; compared with the statistical model checking tool PlasmaLab, HooVer can require much less data to achieve comparable results.
  Regret Bounds and Reinforcement Learning Exploration of EXP-based Algorithms
  EXP-based algorithms are often used for exploration in multi-armed bandi...
  09/20/2020 ∙ by Mengfan Xu, et al. ∙ 9 ∙ share
  Satisficing is a relaxation of maximizing and allows for less risky deci...
  12/23/2015 ∙ by Paul Reverdy, et al. ∙ 0 ∙ share
  In recent years probabilistic model checking has become an important are...
  03/08/2021 ∙ by Alberto Termine, et al. ∙ 0 ∙ share
  Narendra-Shapiro (NS) algorithms are bandit-type algorithms that have be...
  02/17/2015 ∙ by Sébastien Gadat, et al. ∙ 0 ∙ share
  With the availability of big medical image data, the selection of an ade...
  The multi-armed bandit problem is an idealized mathematical model for sequential decision making in unknown random environments and it has been used to study exploration-exploitation trade-offs. In the problem setup, each arm x∈X of the bandit is associated with a cost λx of playing and an unknown reward distribution Mx. In order to maximize the final reward with a given cost budget, the algorithm plays some arm, collects the stochastically generated reward, and decides on the next arm, until the cost budget is exhausted. Starting from the motivation of designing clinical trials in the 1930s [thompson1935theory, thompson1933likelihood, robbins1952some] , there has been major developments in the Bandit theory over the last few decades (see, for example the books [munos:hal-2014, Bubeck:2011, Bubeck12] ). Several different strategies have addressed this problem and strong connections have been drawn with other fields such as online learning.
  In this paper, we explore how Bandit algorithms can be used for model checking of stochastic systems. A requirement R for a stochastic system M usually checks whether the measure of executions of M satisfying certain temporal formulas cross certain thresholds [younes2002probabilistic, GburekB18] . Model checking for such requirements can be solved by calculating the exact measure of the executions that satisfy the relevant subformulas of R [BustanRV04, HermannsWZ08, JansenKOSZ07, Kwiatkowska:book2004] . In this paper, we focus on the alternative statistical model checking (SMC) approach which samples some executions of M and uses hypothesis testing to infer whether the samples provide statistical evidence for the satisfaction (or violation) of R [younes2002probabilistic, Sen:2005:SMC, Younes05] . Execution data is a costly resource111Generating execution data involves running simulations or performing tests., therefore, a number of SMC approaches minimize the expected number of samples needed for verification, for example, using sequential probability ratio tests, Chernoff bound, and Student's t-distribution.
  Several SMC tools have been developed (for example, Ymer [Ymer] , VESTA [VESTA:SenVA05] , MultiVesta [MVesta] , PlasmaLab [Boyer:2013:PFD] , MODES [MODES] , UPPAAL [UPPAAL] , and MRMC [MRMC] ), and they have been used to verify many systems [DDavidLLMPVW11, BaierHHK02, JhaCLLPZ09, KyleHC15, PalaniappanG0HT13, ZulianiPC13, Meseguer:2006:SAD, Martins:2011:SMC]
  . Most SMC algorithms crucially rely on fully stochastic models that never make nondeterministic choices. Although recent progress has been made towards verifying Markov Decision Processes with restricted types of schedulers
  [lassaigne2015approximate, DBLP:conf/tacas/HartmannsH14, henriques2012statistical] , SMC for MDPs remain a challenge problem (see [Agha:2018:SSM] and [Legay:2015] for recent surveys).
  We will focus on stochastic models that are essentially Discrete Time Markov Chains, except that they are initialized from a (possibly very large) set of states. In other words, these are Markov Decision Processes (MDPs) where the adversary gets to initialize222 Finite number nondeterministic action choices can be encoded in the choice of the initial state.. Further, we restrict our attention to safety requirements333All the results in the paper generalize to bounded time properties of the form P≥θ(ψ) where θ is a threshold constant and ψ is a path formula. Generalizing to nested probabilistic operators and unbounded time properties will require further research.
  . That is, we study problems that require maximizing (or minimizing) the probability of hitting certain unsafe states, starting from any initial state. Further, this class of models and requirements capture many practical problems like online monitoring where the initialization has to consider worst case error bounds in state estimation, for example, from sensing and perception. We observe that this optimization of a probability measure over a set of initial choices, can coincide with the multi-armed bandit problem for appropriately defined costs and rewards. By building the connection with the Bandit literature, we not only gain algorithmic ideas, but also new types of theoretical (regret) bounds on the sample efficiency of the algorithms. These bounds rely on the smoothness and the near-optimality dimension of the objective function, and are fundamentally different from the existing performance bounds in the SMC literature.
  Hierarchical optimistic optimization (HOO) [Bubeck:2011] is a bandit algorithm that builds a tree on a search space X by using the so called principle of optimism in the face of uncertainty. It is a black-box optimization method that applies an upper confidence bound (UCB) on a tree search method for finding the optimal points over the uncertain domain. The UCB in the tree search approach takes care of the trade-off between exploiting the most promising parts of the domain and exploring the most uncertain parts of the domain. Multi-fidelity hierarchical optimistic optimization(MFHOO) [sen2019noisy] is a multi-fidelity HOO based method that allows noisy and biased observations from the uncertain domain. The performance of MFHOO is measured by how the regret—the gap between the actual maximum and the computed—scales with the number of samples. A key feature of these algorithms is that they can work with black-box functions and the regret guarantees only rely on certain smoothness parameters and the near-optimality dimension of the problem (see Definition 3.1).
  The standard theoretical assumptions required by off-the-shelf bandit algorithms in order to get performance guarantees do not precisely fit our verification problem, and that in-depth analysis and modification is required to get these guarantees in our setting; In addition, to apply these algorithms several functions need to be judiciously determined, a priori, and are at the heart of how the algorithms will perform. These choices are non-trivial and multi-faceted, and we develop and provide such functions explicitly in the context of our SMC problem in order to demonstrate successful application.
  The key contributions of the paper are as follows.
  First, we show how the MFHOO algorithm, can be used for statistical model checking with provable regret bounds. In the process, we define an appropriate notions of fidelity, bias-functions, and also modify the existing near-optimality dimension required for regret bounds of MFHOO to accommodate the non-smoothness of the typical functions we have to optimize for SMC.
  Second, we have built a new SMC tool called HooVer using MFHOO [sen2019noisy] . We have created a practically inspired [simonepresent, FanQM18] suite of benchmark NiMC models that can be useful for safety analysis of driver assistance features in vehicles for standards such as ISO26262 [iso26262] . Using the benchmarks we have carried out a detailed performance analysis of HooVer and our results suggest that the proposed approach can indeed make use of simulations more effectively than existing SMC approaches. A fair comparison of HooVer with other discrete-state model checkers like Prism [HKNP06] , Storm [dehnert2017storm] , and PlasmaLab [legay2016plasma] is complicated as it relies on a continuous state models. We created discretized models for comparison, and observed that: Compared with exact probabilistic model checking tools like Storm, HooVer is faster, more memory efficient and scales better, and thus it can be used to check models with very large initial state space; Compared with statistical model checking tools like PlasmaLab, HooVer requires much less data to achieve comparable results.
  Finally, to our knowledge, this is the first work connecting statistical model checking with the Bandits theory; specifically, the hierarchical tree search using the principle of optimism in the face of uncertainty. Thus we believe that the exposition of these algorithms (Section 3) engender new applications in verification and synthesis algorithms.
  Consider a Euclidean space X=Rm and let R≥0 denote the non-negative real numbers. For any real-valued function p of X its support is the set supp(p)\coloneqq{x∈X | p(x)≠0}. A discrete probability distribution over X is a function p:X→[0,1] such that supp(p) is countable, and ∑x∈supp(p)p(x)=1. We use P(X)
  to denote the set of discrete probability distributions over
  X. For a finite set S, |S| denotes the cardinality of S.
  2.1 Nondeterministically initialized Markov chains
  A Nondeterministically initialized Markov chains (NiMC) M is defined by a triple (X,Θ,P), where: (i) X=Rmis the state space; (ii) Θ⊆Xis the set of possible initial states; and (iii) P:X→P(X)is the probability transition function.
  That is, from state x∈X, the next state is chosen according to the discrete distribution P(x). The probability of transitioning from state x to state x′∈X is P(x)(x′), which we write more compactly as Px,x′. An execution α of length k for the NiMC M is a sequence of states α={x0,x1,…,xk}, where x0∈Θ, and for each i>0, Pxi−1,xi>0. We denote the set of all length k executions of M starting from x0 as Execsx0(k). The probability of an execution α, given x0, is ∏ki=1Pxi−1,xi=:p(α).
  Given a set U⊆X, we say an execution α hits U if there exists x∈α such that x∈U. We denote the subset of executions starting from x0, of length k, that hit U by Execsx0(k,U). From a given initial state x0∈Θ the probability of M hitting an unsafe state within k steps is given by:
  pk,U(x0)=∑α∈Execsx0(k,U)p(α). (1)
  Note that if x0∈U then Execsx0(k,U)=Execsx0(k) and pk,U(x0)=1. We are interested in finding the worst case probability of hitting unsafe states from any initial state of M. This can be regarded as determining, for each k, the value
  maxx0∈Θpk,U(x0). (2)
  Importantly, we would like to solve this optimization problem without relying on detailed information about the probability transition function P. Further, our solution should not rely on precisely computing pk,U(x0) for a given x0∈Θ, but instead only the use of noisy observations.
  2.2 Example: Single-lane platoon with two speeds (SLplatoon2)
  We present an NiMC of a platoon of m cars on a single lane (SLplatoon2). Variations of this model are used in all our experiments later in Section 5. Each car probabilistically decides to "cruise" or "brake" based on its current gap with the predecessor. These types of models are used for risk analysis of Automatic Emergency Braking (AEB) systems [simonepresent, FanQM18] . The probabilistic parameters of the model are derived from data collected from overhead traffic enforcement cameras on roads. The uncertainty in the initial positions (and gaps) arise from perception inaccuracies, which are modeled as worst-case ranges.
  ⬇ 1if i=1 si←si+vc 3else gapi←si−1−si 5 if gapi≤near si←si+vc w.p. pnear 7 si←si+vb w.p. 1−pnear else 9 si←si+vc w.p. pfar si←si+vb w.p. 1−pfar
  Figure 1: Left: Model of SLplatoon2. Right: Probability of hitting U vs. the initial gap with different execution lengths k. Here, initial s1∈(20,23),s2∈(0,2), vb=1, vc=2, pfar=0.85, pnear=0.15, near=3 and δ=1
  Let si be the position of ith car in the sequence. Initially, si takes a value in an interval on the x-axis such that s1>s2>…>sm. The pseudocode in Figure 1 specifies the probabilistic transition rule that updates the position of all the cars synchronously. Car 1 always moves at a constant breaking speed of vc. The variable gapi is the distance of i to the predecessor i−1, for each i=2,...,m. If gapi is less than the constant threshold near, then i continues to cruise with probability pnear and it brakes with probability 1−pnear. Similarly, if gapi greater than near, then i continues to cruise with probability pfar and brakes with probability 1−pfar.
  It is straightforward to connect the above description to a formal definition of a NiMC. The state space X=Rm. The set of initial states Θ is a hyperrectangle in X (such that s1>s2…>sm). For any state, x∈X the probability transition function is given by the equations in lines 1, 1, 1, 1 and 1. We define the set of unsafe states U={(s1,s2,...,sm)∈X | ∃i∈{2,...,m}, such that gapi≤δ}⊆X for some constant collision threshold δ.
  Given that cars start their motion at any initial state from Θ, the goal is to find the maximum probability of hitting the unsafe set U. For m=2, vb=1, vc=2, pfar=0.85 pnear=0.15, near=3, δ=1, and initial s1∈(20,23),s2∈(0,2), Figure 1 shows estimates of probabilities of hitting the unsafe set from different initial separations between cars. As our intuition suggests, for large enough time horizons the probability of hitting the unsafe set approaches 1 from all initial states, but, for smaller time horizons the maximum probability of unsafety arises when the initial gap is smaller.
  3 Background: Hierarchical Optimistic Optimization
  Multi-Fidelity Hierarchical Optimistic Optimization (MFHOO) [sen2019noisy] is an black-box optimization algorithm from the multi-armed bandits literature [munos:hal-2014, Bubeck:2011, Bubeck12] . The setup is the following: suppose we want to maximize the function f:X→R, which is assumed to have a unique global maximum. Let f∗=supx∈X f(x). MFHOO allows the choice of evaluating f at different fidelities with different costs. This flexibility matters for SMC because it will be beneficial to evaluate the probability of unsafety pk,x0(U) for certain initial states more precisely, for example, with longer number of simulations, while for other initial states a less precise evaluation may be adequate.
  Thus, MFHOO has access to a biased function fz(x) that depends on fidelity parameter z∈[0,1]. Setting z=0 gives the lowest fidelity (and lowest cost) and z=1 corresponds to full fidelity (and highest cost). At full fidelity, f1(x)=f(x), and the evaluation is unbiased. More generally, |fz(x)−f(x)|≤ζ(z) and evaluating fz(x) costs λ(z), where the functions ζ,λ:[0,1]→R>0 are respectively, non-increasing and non-decreasing, and called the bias and the cost functions [sen2019noisy] .
  A bandit algorithm chooses a sequence of sample points (arms) x1,x2,…∈X, evaluates them at fidelities z1,z2,…, and receives the corresponding sequence of noisy observations (rewards) y1,y2,…. We assume that each yj is drawn from a unknown distribution Mzj,xj satisfying fzj(xj)=∫udMzj,xj(u)
  . Further the distribution has a sub-Gaussian component, with variance
  σ2, which captures uncertainty in the observations. The algorithm actively chooses xj+1 based on past choices x1,…,xj and observations y1,…,yj. When the budget Λ is exhausted, the algorithm decides the optimal point ¯xn(Λ)∈X and the optimal value f(¯xn(Λ)) with the aim of minimizing regret, which is defined as
  R(Λ)=f∗−f(¯xn(Λ)).
  The MFHOO algorithm (Algorithm 1) for selecting xj+1 estimates f∗ by building a binary tree in which each height in the tree represents a partition of the state space X. The algorithm maintains estimates of an upper-bound on f for each partition subset, and uses the principle of optimism for choosing the next sample xj+1. That is, it chooses the samples in those partitions where the estimated upper-bounds are the highest.
  Each node in the constructed tree is labeled by a pair of integers (h,i), where h is the height of the node in the tree, and i satisfying 0≤i≤2h is its position within height level h. The root is labeled (0,1), and each node (h,i) can have two children (h+1,2i−1) and (h+1,2i). Node (h,i) is associated with subset a of X, denoted by Ph,i, where Ph,i=Ph+1,2i−1∪Ph+1,2i, and for each h these disjoint subsets satisfy ∪2hi=1Ph,i=X. Therefore, larger values of h represent finer partitions of X.
  For each node (h,i) in the tree, the algorithm maintains the following information: (i) counth,i: the number of times the node is visited; (ii) ^fh,i: the empirical mean of observations over points visited in Ph,i; (iii) Uh,i: an initial estimate of the maximum of f over Ph,i; and (iv) Bh,i: a tighter and optimistic upper bound on the maximum of f over Ph,i.
  The algorithm proceeds as follows. The tree starts out with a single node, the root (0,1), initializes the B-values of its two children B1,1 and B1,2 to +∞, and initializes the cost C to 0. At a high-level, in each iteration of the while loop (line 3), the algorithm adds a new node (hnew,inew) in the tree and updates all of the above quantities for several nodes in tree. In more detail, first a path from the root to a leaf is found by traversing the child with the higher B-value (with ties broken arbitrarily). Let the child with the higher B-value of the traversed leaf be (hnew,inew) (line 4). An arbitrary point x in the partition of this node Phnew,inew is chosen (line 5). Then, this point is evaluated at fidelity zhnew=ζ−1(νρhnew) and a reward y is received (line 6). Next, tree is extended by inserting (hnew,inew) in the tree (line 8) and for all the nodes (h,i) in path including (hnew,inew), the counth,i and the empirical mean ^fh,i are updated (line 9). Finally, in line 13, for all nodes (h,i) in tree, Uh,i and Bh,i are updated using the smoothness parameters ν1>0 and ρ∈(0,1) which will discussed later in Section 3.1 and the parameter σ. Once the sampling budget Λ is exhausted, a leaf with maximum B-value is returned by the Algorithm 1 [sen2019noisy] .
  1:input: Budget: Λ, parameter: σ, bias ζ(.), cost λ(.), smoothness params ν>0 and ρ∈(0,1)
  2:tree={(0,1)}, B1,1=B1,2=∞, C=0
  3:while C≤Λ do
  4: (path,(hnew,inew))←Traverse(tree)
  5: choose x∈Phnew,inew
  6: query x at fidelity zhnew=ζ−1(νρhnew) and get observation y
  7: C←C+λ(zhnew)
  8: tree.Insert((hnew,inew))
  9: for all (h,i)∈path do
  10: counth,i←counth,i+1
  11: ^fh,i←(1−1/counth,i)^fh,i+y/counth,i
  12: Bhnew+1,2inew−1←+∞, Bhnew+1,2inew←+∞
  13: for all (h,i)∈tree do from leaves up to root:
  14: Uh,i←^fh,i+sqrt(2σ2lnn/counth,i)+νρh+ζ(zh)
  15: Bh,i←min{Uh,i,max{Bh+1,2i+1,Bh+1,2i}}
  16:return argmax(h,i)∈tree Bh,i
  Algorithm 1 Multi-Fidelity Hierarchical Optimistic Optimization [sen2019noisy]
  3.1 Regret bounds for MFHOO
  In this section, we summarize the assumptions and results from [sen2019noisy] . In order to analyze the regret bounds for the MFHOO algorithm, the following assumption on the smoothness of f(x) [sen2019noisy] is made.
  There exist ν>0 and ρ∈(0,1) such that for all Ph,i satisfying f∗−supx∈Ph,i f(x)≤cνρh (for a constant c≥0), we have that f(x)≥f∗−max{2c,c+1}νρh , for all x∈Ph,i.
  This assumption connects the function f to the partitioning rule of the binary tree. We now define the concept of near-optimality dimension which plays an important role in the analysis of black-box optimization algorithms [Bubeck:2011, sen2019noisy, grill2015black, valko2013stochastic] . It measures the dimension of sets that are close to optimal. The regret bound for MFHOO uses this notion. First, given a partitioning scheme Ph,i over X, we define Nh(ϵ) as the number of ϵ-near optimal partitions, that is, the number of partitions Ph,i such that satisfies supx∈Ph,if(x)≥f(x∗)−ϵ.
  The near-optimality dimension of f with respect to parameters (ν,ρ) is: d(ν,ρ)=inf{d′∈R>0:∃B>0, s.t. ∀h≥0, Nh(2νρh)≤Bρ−d′h}.
  With Assumption 1, the regret bound for MFHOO is proved in [sen2019noisy] .
  If Algorithm 1 runs with parameters ν and ρ that satisfy Assumption 1 and a cost budget of Λ, then the simple regret is bounded
  R(Λ)=O((Blogn(Λ)n(Λ))1d(ν,ρ)+2),
  where n(Λ)=max{n:∑nh=1λ(zh)≤Λ}. Here, zh=ζ−1(νρh).
  According to the Theorem 1, regret is minimized if the near-optimality dimension d(ν,ρ) is minimized. If the smoothness parameters that minimize the near-optimlaity dimension d(⋅,⋅) are known, then MFHOO achieves the minimum regret of Theorem 1.
  4 Statistical Model Checking with Optimistic Optimization
  We aim to solve the statistical model checking problem of maximizing pk,U(x) of Equation (2) for a given NiMC M and a time horizon k, using MFHOO. In order to apply the MFHOO algorithm, one has to make several critical choices regarding the objective function, the budget, the cost, the parameters for fidelity and smoothness, and the multi-fidelity bias function. In this section we discuss the rationale behind our choices.
  4.1 Objective function, budget, cost, and fidelity
  Fidelity parameter z.
  Consider a NiMC M=(X,Θ,P) with the unsafe set U⊆X. We have to maximize pkmax,U(x) over all initial states x∈Θ, and for a long time horizon kmax. Given x∈Θ, the fidelity of evaluating pkmax,U(x) will depend on the actual length of the simulations drawn for creating the observation y for the state x. Suppose we fix kmin as the shortest simulation to be used. We define the fidelity of an observation (or evaluation) with simulations of length k∈[kmin,kmax] as z=(k−kmin)/(kmin−kmax).
  Objective function f and observations.
  A natural choice for the objective function would be to define fz=1(x):=pkmax,U(x), for any initial state x∈Θ. Computing this probability, however, is infeasible when the probability transition function PM is unknown. Even if PM is known, calculating pk,U(x) involves summing over many executions (as in (1)). Instead, we take advantage of the fact that MFHOO can work with noisy observations. For any initial x∈Θ, and execution α∈Execsx(k) we define the observation:
  Y=1 if α∈Execsx(k,U), and =0 otherwise. (3)
  Recall that for a fixed initial states x, M is a Markov chain and defines a probability distribution over the set of executions Execsx(k) as given by Equation (1). Thus, given an initial state x, Y=1 with probability pk,U(x), and Y=0 with probability 1−pk,U(x). That is, Y
  pk,U(x) at fidelity z. In MFHOO, once an initial state x∈Phnew,inew is chosen (line 5), we simulate M upto k steps several times starting from x and calculate the empirical mean of Y, which serves as the noisy observation y at fidelity z.
  Cost λ(z) and budget Λ.
  In our setup the cost function λ(z) for any z is the computational time required to simulate an execution α∈Execsx(k) for any x∈Θ, where k is the execution length corresponding to the fidelity z. The budget Λ is a computational cost budget.
  The next proposition states that with these choices, our goal to maximize pkmax,U(x) over Θ can be achieved using MFHOO.
  Given smoothness parameters ρ and ν satisfying Assumption 1 and budget Λ, suppose Algorithm 1 returns ¯xΛ∈Θ. Then, R(Λ)=pkmax,U(x∗)−pkmax,U(¯xΛ), where R(Λ) bound is given by Theorem 1.
  Thus, the point ¯xΛ∈Θ returned by Algorithm 1 is such that the gap between its probability pk,U(¯x) of hitting U and the true maximum probability pk,U(x∗) is bounded by Theorem 1 in terms of the available budget Λ.
  4.2 Multi-fidelity bias function
  Recall that the bias function ζ(z) gives an upper bound |f(x)−fz(x)|≤ζ(z), over all x∈Θ and for any fidelity z∈[0,1]. We will derive a bias function satisfying |pkmax,U(x0)−pk,U(x0)|≤ζ(z). Of course, a guarantee like this is only possible for known models. Therefore, for this section we will assume that the NiMC M is known. We also assume X is finite and U is the absorbing set for M (i.e., all other states are transient; x∈U if and only if Pxx=1).
  Fixing an initial state x, M is a reducible Markov chain. Let q be the number of transient states and u=|U|. Then, the probability transition function can be represented by the (q+u)×(q+u) matrix P:
  P=(QR0I) and Pt=(Qt∑t−1k=0QkR0I). (4)
  Here Q is an q×q matrix, I is an u×uidentity matrix, 0\/ is an u×qzero matrix, and R is a nonzero q×u matrix. The following standard result gives the absorption probabilities in terms of Q.
  Proposition 2 ( [grinstead2012introduction] )
  Suppose B=NR is a q×u matrix, where N=∑∞j=0Qj=(I−Q)−1 and R is as defined in (4). Then Bij is the probability, that starting from state si, the chain M is absorbed in sj∈U.
  We now state and prove the theorem that defines a multi-fidelity bias function.
  For any x∈Θ and any time horizon k<kmax,
  |pkmax,U(x)−pk,U(x)|≤κ∞(E)(kmax−k)(kmax−(kmax−1)λmax)(k−(k−1)λmax),
  is the maximum eigenvalue,
  is the matrix of eigenvectors of the matrix
  Q, and κ∞(E) is the (∞-norm) condition number of E.
  This theorem can be proved using the spectral decomposition of matrix Q and ∞-norm bounds (for detailed analysis please refer to Appendix B).
  Using the fidelity parameter z=(kmax−k)/(kmax−kmin) in the upper bound given by Theorem 2, the bias function can be rewritten as:
  ζ(z) =g1(1−z)g2z+g3,where g1=κ∞(E)(kmax−kmin) (5)
  g2 =((1−λmax)kmax−λmax)(kmax−kmin)(1−λmax)
  g3 =((1−λmax)kmax−λmax)((1−λmax)kmin−λmax).
  This bias function upper bounds the gap between the probability of hitting the unsafe set for a point in the initial set for different time horizons. Thus, it can be used in the analysis of the theoretical regret bound of the Algorithm 1 given in Theorem 1. If the transition matrix P is known, then the bias function can be used in updating the U-values of nodes in the Algorithm 1. Note, this bias function depends on parameters κ∞(E) and λmax. In the case where the full transition model of the NiMC is unknown, but there is access to κ∞(E) and λmax, this function can be utilized in the algorithm. More generally, in problems with unknown P and no access to the parameters κ∞(E) and λmax, we consider a linear parameterized bias function ζ(z)=b(1−z) with an unknown parameter b, and adaptively estimate the parameter b [sen2019noisy] .
  4.3 Non-smoothness and non-uniqueness of hitting probabilities
  In this section, we start with the observation that in general the function pk,U(x) over x∈Θ is not a continuous function and has infinitely many maxima. Thus this function does not satisfy the Assumption 1 and the assumption of finite number of maxima, that are required for the regret bounds of Theorem 1. However, a modified definition of near-optimality dimension we show that the bounds given by Theorem 1 will hold. Finally, we will compare the theoretical regret with the actual regret achieved by MFHOO.
  Consider a discrete time linear system with state space X=Rm, set of initial states Θ⊆X and set of unsafe states given by U:={x∈X | cTx≤a}
  c∈Rm and a∈R. The system evolves according to:
  xt+1=Axt+BWt, (6)
  where xt∈X is the state vector at time t. Wt∈Rr is a discrete random vector with probability mass distribution pW at sequence t, and A and B are m×m and m×r matrices, respectively.
  Consider two initial states x0,x′0∈Θ, such that x′0=x0+ϵ, for some small ϵ. Writing out the probabilities of hitting unsafe states of this system explicitly for the two initial states and finding the gap between these probabilities we can observe that, for small enough ϵ, the gap between the probability is zero444 The full analysis is given in Appendix D.. By increasing the ϵ, the gap can take a nonzero value. This means that the probability of hitting the unsafe set of the system (13) is not a continuous function in terms of the initial states and has infinitely many maxima. This is because random variable Wt is a discrete-type random variable. SLplatoon2 is a concrete example of this kind.
  Assumption 1 asserts that there exist smoothness parameters such that for all h≥0 in the tree, and all the partitions that are near optimal, the gap between f∗ and the value of the function over those partitions should be bounded. From the above discussion we see that this may not hold for all depths h≥0. In addition, the regret bounds in Theorem 1 are for function with finite number of maxima. If there are infinitely many maxima, then for given smoothness parameters the number of 2νρh-near-optimal partitions Nh(2νρh) would not be bounded for all h≥0. However, we observe that any instance of MFHOO runs on a finite budget and the final constructed tree has a maximum height hmax. For this hmax, there exist smoothness parameters ν and ρ that satisfy the Assumption 1, and we can modify the near-optimality dimension in the Definition 3.1 as:
  hmax-bounded near-optimality dimension of f with respect to (ν,ρ) is: dm(ν,ρ)=inf{d′∈R>0:∃B>0, ∀h∈[0,hmax],Nh(2νρh)≤Bρ−d′h}.
  With this modified definition, there exists a B satisfying the above condition and the corresponding dm can be used to recover the regret bound of the Theorem 1.
  Regret in theory and in practice.
  Given a budget Λ, let ¯xΛ∈Θ be the point that the Algorithm 1 returns after the budget is exhausted. Then the regret R(Λ) would be pkmax,U(x∗)−pkmax,U(¯xΛ). We compare the theoretical regret bound and the actual regret for the SLplatoon2 model. Consider an instance of SLplatoon2 with m=2, initial states s1∈(50,70) and s2∈(0,20). Intuitively, the partition corresponding to s1∈(50,51) and s2∈(19,20) would maximize the probability of hitting U. In Figure 2 (top-left), the mean of actual regret obtained by running MFHOO for various smoothness parameters is presented. As the budget increases, as expected, the actual regret decreases monotonically and approximately approaches to 0.
  For budget Λ=30, the partitioning stops at hmax=15. We can derive the number of 2νρh-near-optimal partitions for any ν and ρ. Setting ν=0.08 ρ=0.77 (actual values used to generate regret results), we see that for h∈[9,15], the maximum number of 2νρh-near-optimal partitions Nh(2νρh)=⌈2h/400⌉. This is because the total number of partitions at depth h is equal to 2h, and for h∈[9,15], the 2νρh-near-optimal partitions belong to the area corresponding to s1∈(50,51) and s2∈(19,20), whose area is 1/400. According to Definition 4.1, for different values of B, different hmax-bounded near-optimality dimensions dm can be obtained. We are looking for the pair of (dm,B) values that minimize the theoretical regret bound of Theorem 1. Figure 2 (top-right) represents the number of 2νρh-near-optimal partitions for h∈[9,15] that are upper bounded by Bρ−dmh for different values of dm and their corresponding B. Figure 2 (bottom) also, represent the theoretical upper bound vs. dm and their corresponding B for different values of smoothness parameters used to generate the actual regret. As it is seen, the best theoretical regret bounds that can be achieved is approximately 0.08.
  Figure 2: SLplatoon2 with the same parameters as in Section 2.2 and initial states s1∈(50,70) and s2∈(0,20). Top-left: Actual regret of running MFHOO averaged for various smoothness parameters. Top-right: Number of 2νρh partitions for various dm values and their corresponding B with ν=0.08 and ρ=0.77. Bottom-left: Plot of g(dm) vs. dm for various smoothness parameters, where R(Λ)=O(g(dm)) as in the Theorem 1. Bottom-right: Plot of g(B) vs. B for various smoothness parameters, where R(Λ)=O(g(B)) as in the Theorem 1.
  5 HooVer tool and experimental evaluation
  We have implemented a prototype tool called HooVer which uses MFHOO for solving the SMC problem of (2). We compare the performance of HooVer with that of Prism [HKNP06] , Storm [dehnert2017storm] , and PlasmaLab [legay2016plasma] on several benchmarks we have created.
  We have created several NiMC models for evaluation of probabilistic and statistical model checking tools. The benchmarks are variants of SLplatoon2 (Section 2.2). The complete models are described in Appendix C. The executable models are available from [SMCOO-online:2020] .
  SLplatoon3 models a sequence of m cars on a single lane. At each step, each car can choose to move with one of three speeds: vbrake, vcruise, and vspeedup. The first vehicle moves at a constant speed. For all the others, the speed is chosen probabilistically, according to distributions that depend on their distance to the preceding vehicle. For example, if the longitudinal distance si−1−si is less than a threshold th_close then the speed for vehicle i is chosen according to a probability distribution pclose, and so on. The state variables, the initial states, and the unsafe set are defined in the same way as in SLplatoon2.
  MLplatoon models m cars on ℓ lanes. At every step, each car probabilistically chooses to either move forward like a vehicle in SLplatoon3, or it changes to its left or right lane. These actions are chosen probabilistically, according to probability distributions that depend on their distances to the vehicles on its current lane, as well as left and right lanes. In addition to the longitudinal position si, each vehicle i has a second state variable yi∈{1,…,ℓ} which keeps track of its current lane. The initial value of each si is chosen as in the case of SLplatoon2, and yi is set to 1. The unsafe set is defined based on the distance to the preceding car on the same lane.
  5.2 HooVer implementation and metrics
  Our implementation of the HooVer tool uses the MFHOO implementation presented in [sen2019noisy] to solve the model checking problem of Equation (2). It works in two stages: First, it uses MFHOO to find the best partition Ph,i and a putative "best" (most unsafe) initial point x∈Ph,i with the maximum estimate for the probability of hitting the unsafe set U. In the second stage, HooVer uses additional simulations to do a Monte Carlo estimation of the probability pk,U(x). In the experiments reported below, a constant number of 26K simulations are used in all experiments in the second stage.
  To achieve the theoretically optimal performance, MFHOO requires the smoothness parameters ρ and ν which are unknown for our benchmarks. To circumvent this HooVer chooses several parameter configurations (3 sets in our experiments), runs an instance of MFHOO for each, and returns the result with the highest hitting probability. For each instance of MFHOO, we set a time budget which is the maximum time allowed to be consumed by the simulator.
  We report the regret of HooVer. In order to calculate the regret, first we have to calculate the actual maximum hitting probability for each benchmark. This is computed using Prism [HKNP06] which uses numerical and symbolic analysis. The regret is the difference between the exact probability and the estimated probability. Then, we report the running time and memory usage. The memory usage is measured by calculating the total size of the Python object which contains the constructed tree and all other data of MFHOO. We also report the number of queries for each method, which is total number of simulations used in both stage 1 and 2. All the experiments are conducted on a Linux workstation with two Xeon Silver 4110 CPUs and 32 GB RAM.
  Table 1 shows the running time, the memory usage, the number of queries, and the resulting actual regret for HooVer using MFHOO as well as HooVer(1) using HOO. On every benchmark, HooVer gives low regrets. With the same simulation budget, HooVer devotes longer simulations in the interesting parts Θ, as a consequence, it is usually faster than HooVer(1) as shown in Figure 3.
  SLplatoon2 m=2 pkmax,U(x∗)=0.8800 Storm 68.88 1019 NA 0
  PlasmaLab 37.73 NA 749711 0.0017 ± 0.0020
  HooVer(1) 59.37 9.13 31381 0.0011 ± 0.0025
  HooVer 24.03 6.23 29415 0.0002 ± 0.0020
  SLplatoon3 m=2 pkmax,U(x∗)=0.8727 Storm 89.10 1974 NA 0
  PlasmaLab 22.61 NA 749711 0.0022 ± 0.0021
  HooVer(1) 44.52 7.30 30315 0.0038 ± 0.0113
  HooVer 26.46 4.53 28520 0.0010 ± 0.0018
  MLplatoon m=2, ℓ=2 pkmax,U(x∗)=0.5918 Storm NA OOM NA NA
  PlasmaLab 49.14 NA 749711 0.0032 ± 0.0026
  HooVer(1) 110.64 9.37 31555 0.0016 ± 0.0025
  HooVer 76.35 5.27 28955 0.0007 ± 0.0043
  Table 1: Comparison of HooVer, Storm and PlasmaLab. HooVer(1) corresponds to using HOO (i.e. the full fidelity algorithm) in stage 1 of HooVer. For regret of HooVer and PlasmaLab, we run the experiments for 10 times and report the "mean ± std". For all the experiments, |¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\Uptheta|=4096.
  5.3 Comparison with other model checkers
  We compare the performance of HooVer with other model checkers. Prism [HKNP06] and Storm [dehnert2017storm] are leading probabilistic model checkers for Markov chains and MDPs and compute exact probability of reaching the unsafe states. As Storm has the same functionality as Prism and we found it to be much more efficient than Prism in all our experiments, we only compare HooVer with Storm. PlasmaLab [legay2016plasma] is one of the few statistical model checkers that can handle MDPs. For probability estimation problems, PlasmaLab uses smart sampling algorithm [d2015smart] to efficiently assign the simulation budgets to each scheduler and then estimates the probability for the putative "best" scheduler.
  Discretizing and scaling the benchmarks.
  The theory for MFHOO is based on a continuous state spaces, however, most model checking tools, including the ones mentioned above, are designed for discrete state space models and they do not support floats as state variables. Therefore, a direct comparison of the approaches on the same model is infeasible, and we created equivalent, discretized (quantized) versions of the benchmarks.
  In HooVer, the algorithm keeps partitioning Θ hierarchically and stops at a depth h, which can be considered as searching over all the 2h partitions at depth h. In order to make a fair comparison, we make sure that the discrete version of the benchmark has 2h initial states, i.e. |¯¯¯¯Θ|=2h where ¯¯¯¯Θ is the discretized initial state space.
  Before stating the discretizing process, we give two key observations of the benchmarks. First, considering the nature of our benchmarks, it is obvious that if we scale the velocities, distance thresholds and initial distances by a constant factor, the probability of hitting unsafe set doesn't change. Second, taking SLplatoon2 as an example, we set all the constants (velocities and distance thresholds) in the model as integers, which leads to the function pkmax,U(x) shown in Figure 1. It's clear that for any state x∈X, there exists an integer state ¯¯¯x such that pkmax,U(x)=pkmax,U(¯¯¯x). If we make sure that for each integer interval in the original continuous space, the discretized space has an value in that interval, then the maximum probability doesn't change. With these observations, we discretize the benchmark as follows. First, we sample 2h points uniformly in the continuous initial state space. Due to the second observation mentioned above, if 2h is larger than the number of integer points in the original space, the maximum probability doesn't change. Then, we find an integer scaling factor c such that by multiplying c all the 2h points become integer points. Other constants in the model are also multiplied by c. According to the first observation, scaling with a constant doesn't change the probability. Thus, we now have a model where all state variables are integers and it has the same maximum probability as the original model. Then we evaluate Storm and PlasmaLab on this new model.
  All of the model checkers mentioned above support the Prism [HKNP06] language. Thus, we implement each benchmark in Prism language, and then we check the equivalence between the Prism implementation and the Python implementation by calculating and comparing the probability pkmax,U(x) at all integer points x.
  For Storm, we report the running time and the memory usage. These metrics are directly measured by the software itself. For PlasmaLab, we report the running time. We do not report the memory usage of PlasmaLab because it doesn't provide an interface for that and it is hard to track the actual memory used by the algorithm inside the JAVA virtual machine. We also report the regret for PlasmaLab. We use the term "regret" here just for simplicity, which also refers to the difference between the estimated probability and the exact probability.
  The results of HooVer, PlasmaLab and Storm are summarized in Table 1. We show in Figure 3 how the performance of different methods changes as the |¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\Uptheta| increases. As the size of the initial state space increases, the memory usage of Storm grows quickly, which limits its application on large models. In contrast, HooVer and PlasmaLab scale well. The running time of PlasmaLab is determined by the parameters of the smart sampling [d2015smart] algorithm. We use the same parameters regardless of |¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\Uptheta|, and thus the running time of PlasmaLab is almost constant.
  As shown in Table 1, Storm fails to check the MLplatoon model due to the memory limitation. Compared with PlasmaLab, HooVer requires more running time. However, we note that PlasmaLab is a considerably more mature tool than HooVer. As shown in Table 1, compared to PlasmaLab, HooVer requires much fewer queries to reach comparable regrets, which attests to the sample efficiency of our proposed method.
  Figure 3: Left: Comparison of HooVer and HooVer(1) on SLplatoon3. (Here, we only consider the number of queries in stage 1, because that for stage 2 is just a constant.) Middle: Performance of different methods as the |¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\Uptheta| increases for SLplatoon3. Right: Performance of different methods as the |¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\Uptheta| increases for MLplatoon.
  In this paper, we formulated the statistical model checking problem for a special type of MDP models as a multi-armed bandit problem and showed how a hierarchical optimistic optimization algorithm from [sen2019noisy] can be used to solve it. In the process, we modified the existing definition of near-optimality dimension to accommodate the non-smoothness of the typical functions we have to optimize and recovered the regret bounds of the algorithm. We created several benchmarks, developed a SMC tool HooVer, and experimentally established the sample efficiency and scalability of the method.
  In order to enlarge the area of application of this method, it is necessary to study more general temporal properties for NiMC beyond bounded time safety. In this paper, we focused on a very special type of MDP models with nondeterminism only in the initial set. The results suggest that it will be interesting to study black-box optimization algorithms for model checking more general MDP models.
  Appendix A Multi-Fidelity Parallel Optimistic Optimization (MFPOO)
  The optimal smoothness parameters are generally not known. Algorithm 2 [sen2019noisy] adaptively searches for the optimal smoothness parameters by releasing several MFHOO instances with various ν and ρ values.
  1:input: νmax and ρmax), Λ(.), λ(.)
  2:Let N=(1/2)Dmaxlog(Λ/logΛ), where Dmax=log2/log(1/ρmax)
  4: Spawn MFHOO with parameters (νmax,ρi=ρN/(N−i−1)max) with budget (Λ−Nλ(1))/N
  5:Let ^xΛ,i be the point returned by the ith MFHOO instance for i∈{0,..,N−1}. Evaluate all {^xΛ,i}i at z=1.
  6:return The point ^xΛ=^xΛ,i∗ where i∗=argmaxi f(xΛ,i)+ϵi.
  Algorithm 2 Multi-Fidelity Parallel Optimistic Optimization (MFPOO)
  The following theorem gives the upper bound for the algorithm 2.
  If Algorithm 2 is run with parameters (νmax≥ν∗), (ρmax≥ρ∗) and given a total cost budget Λ, then the simple regret of at least one of the MFHOO instances spawned by Algorithm 2 is bounded a follows
  R(Λ)=O((νmaxν∗)Dmax(Blogn(Λ/logΛ)n(Λ/logΛ))1d(ν∗,ρ∗)+2).
  It is noted that if the Algorithm 1 is run at full fidelity case z=1, then it will be equivalent to a budgeted HOO algorithm. At this case the simple regret bound will be:
  R(Λ)=O((log(Λ/λ(1))Λ/λ(1))1d(ν,ρ)+2). (7)
  (Theorem 2) We define a new Markov chain M(k) that starts from x (like M) and stops by going entering a new absorbing state called Stop, in k steps, in expectation. That is, from any transient state si, the one step transition probability of hitting Stop is 1k. For any other transient state sj, the one step probability of transitioning from si to sj is (1−1k)Pij. Therefore, the probability transition matrix Pk for the new Markov chain M(k) is as follows.
  Pk=⎛⎜ ⎜⎝QkRk1k% 10I0001⎞⎟ ⎟⎠=⎛⎜ ⎜⎝(1−1k)Q(1−1k)R1k10I0001⎞⎟ ⎟⎠,
  where 1 is a q×1 vector whose elements are 1, and the Stop state is the last state. Pk is a (q+u+1)×(q+u+1)