Sound Value Iteration

Computing reachability probabilities is at the heart of probabilistic model checking. All model checkers compute these probabilities in an iterative fashion using value iteration. This technique approximates a fixed point from below by determining reachability probabilities for an increasing number of steps. To avoid results that are significantly off, variants have recently been proposed that converge from both below and above. These procedures require starting values for both sides. We present an alternative that does not require the a priori computation of starting vectors and that converges faster on many benchmarks. The crux of our technique is to give tight and safe bounds - whose computation is cheap - on the reachability probabilities. Lifting this technique to expected rewards is trivial for both Markov chains and MDPs. Experimental results on a large set of benchmarks show its scalability and efficiency.


Introduction
Markov decision processes (MDPs) [1,2] have their roots in operations research and stochastic control theory. They are frequently used for stochastic and dynamic optimization problems and are widely applicable in, e.g., stochastic scheduling and robotics. MDPs are also a natural model in randomized distributed computing where coin flips by the individual processes are mixed with nondeterminism arising from interleaving the processes' behaviors. The central problem for MDPs is to find a policy that determines what action to take in the light of what is known about the system at the time of choice. The typical aim is to optimize a given objective, such as minimizing the expected cost until a given number of repairs, maximizing the probability of being operational for 1,000 steps, or minimizing the probability to reach a "bad" state.
Probabilistic model checking [3,4] provides a scalable alternative to tackle these MDP problems, see the recent surveys [5,6]. The central computational issue in MDP model checking is to solve a system of linear inequalities. In absence of non-determinism -the MDP being a Markov Chain (MC) -a linear equation system is obtained. After appropriate pre-computations, such as determining the states for which no policy exists that eventually reaches the goal state, the (in)equation system has a unique solution that coincides with the extremal value that is sought for. Possible solution techniques to compute such solutions include policy iteration, linear programming, and value iteration. Modern probabilistic model checkers such as PRISM [7] and Storm [8] use value iteration by default. This approximates a fixed point from below by determining the probabilities to reach a target state within k steps in the k-th iteration. The iteration is typically stopped if the difference between the value vectors of two successive (or vectors that are further apart) is below the desired accuracy ε.
This procedure however can provide results that are significantly off, as the iteration is stopped prematurely, e.g., since the probability mass in the MDP only changes slightly in a series of computational steps due to a "slow" movement. This problem is not new; similar problems, e.g., occur in iterative approaches to compute long-run averages [9] and transient measures [10] and pop up in statistical model checking to decide when to stop simulating for unbounded reachability properties [11]. As recently was shown, this phenomenon does not only occur for hypothetical cases but affects practical benchmarks of MDP model checking too [12]. To remedy this, Haddad and Monmege [13] proposed to iteratively approximate the (unique) fixed point from both below and above; a natural termination criterion is to halt the computation once the two approximations differ less than 2·ε. This scheme requires two starting vectors, one for each approximation. For reachability probabilities, the conservative values zero and one can be used. For expected rewards, it is non-trivial to find an appropriate upper bound -how to "guess" an adequate upper bound to the expected reward to reach a goal state? Baier et al. [12] recently provided an algorithm to solve this issue. This paper takes an alternative perspective to obtaining a sound variant of value iteration. Our approach does not require the a priori computation of starting vectors and converges faster on many benchmarks. The crux of our technique is to give tight and safe bounds -whose computation is cheap and that are obtained during the course of value iteration -on the reachability probabilities. The approach is simple and can be lifted straightforwardly to expected rewards. The central idea is to split the desired probability for reaching a target state into the sum of (i) the probability for reaching a target state within k steps and (ii) the probability for reaching a target state only after k steps.
We obtain (i) via k iterations of (standard) value iteration. A second instance of value iteration computes the probability that a target state is still reachable after k steps. We show that from this information safe lower and upper bounds for (ii) can be derived. We illustrate that the same idea can be applied to expected rewards, topological value iteration [14], and Gauss-Seidel value iteration. We also discuss in detail its extension to MDPs and provide extensive experimental evaluation using our implementation in the model checker Storm [8]. Our experiments show that on many practical benchmarks we need significantly fewer iterations, yielding a speed-up of about 20% on average. More importantly though, is the conceptual simplicity of our approach.

Preliminaries
For a finite set S and vector x ∈ R |S| , let x[s] ∈ R denote the entry of x that corresponds to s ∈ S. Let S ′ ⊆ S and a ∈ R. We write x[S ′ ] = a to denote that For a function f : R |S| → R |S| and k ≥ 0 we write f k for the function obtained by applying f k times, i.e., f 0 (

Probabilistic Models and Measures
We briefly present probabilistic models and their properties. More details can be found in, e.g., [15].
Definition 1 (Probabilistic Models). A Markov Decision Process (MDP) is a tuple M = (S, Act, P, s I , ρ), where -S is a finite set of states, Act is a finite set of actions, s I is the initial state, -P : S × Act × S → [0, 1] is a transition probability function satisfying s ′ ∈S P(s, α, s ′ ) ∈ {0, 1} for all s ∈ S, α ∈ Act, and ρ : S × Act → R is a reward function. M is a Markov Chain (MC) if |Act| = 1. Example 1. Fig. 1 shows an example MC and an example MDP.
We often simplify notations for MCs by omitting the (unique) action. For an MDP M = (S, Act, P, s I , ρ), the set of enabled actions of state s ∈ S is given by Act(s) = {α ∈ Act | s ′ ∈S P(s, α, s ′ ) = 1}. We assume that Act(s) = ∅ for each s ∈ S. Intuitively, upon performing action α at state s reward ρ(s, α) is collected and with probability P(s, α, s ′ ) we move to s ′ ∈ S. Notice that rewards can be positive or negative.
A state s ∈ S is called absorbing if P(s, α, s) = 1 for every α ∈ Act(s). A path of M is an infinite alternating sequence π = s 0 α 0 s 1 α 1 . . . where s i ∈ S, α i ∈ Act(s i ), and P(s i , α i , s i+1 ) > 0 for all i ≥ 0. The set of paths of M is denoted by Paths M . The set of paths that start at s ∈ S is given by Paths M,s . A finite path π = s 0 α 0 . . . α n−1 s n is a finite prefix of a path ending with last (π) = s n ∈ S. |π| = n is the length ofπ, Paths M fin is the set of finite paths of M, and Paths M,s fin is the set of finite paths that start at state s ∈ S. We consider LTL-like notations for sets of paths. For k ∈ N ∪ {∞} and G, H ⊆ S let H U ≤k G = {s 0 α 0 s 1 · · · ∈ Paths M,sI | s 0 , . . . , s j−1 ∈ H, s j ∈ G for some j ≤ k} denote the set of paths that, starting from the initial state s I , only visit states in H until after at most k steps a state in G is reached. Sets H U >k G and H U =k G are defined similarly. We use the shorthands ♦ ≤k G := S U ≤k G, ♦G := ♦ ≤∞ G, and ≤k G := Paths M,sI \ ♦ ≤k (S \ G).
A (deterministic) scheduler for M is a function σ : Paths M fin → Act such that σ(π) ∈ Act(last (π)) for allπ ∈ Paths M fin . The set of (deterministic) schedulers for M is S M . σ ∈ S M is called positional if σ(π) only depends on the last state ofπ, i.e., for allπ,π ′ ∈ Paths M fin we have last (π) = last (π ′ ) implies σ(π) = σ(π ′ ). For MDP M and scheduler σ ∈ S M the probability measure over finite paths is given by Pr M,σ fin : The probability measure Pr M,σ over measurable sets of infinite paths is obtained via a standard cylinder set construction [15]. For k ∈ N∪{∞}, the function ≤k G : ♦G → R yields the k-bounded reachability reward of a path π = s 0 α 0 s 1 · · · ∈ ♦G. We set ≤k G(π) = We write Pr M,σ s and E M,σ s for the probability measure and expectation obtained by changing the initial state of M to s ∈ S. If M is a Markov chain, there is only a single scheduler. In this case we may omit the superscript σ from Pr M,σ and E M,σ . We also omit the superscript M if it is clear from the context. The maximal reachability probability of M and G is given by Pr max (♦G) = max σ∈S M Pr σ (♦G). There is a a positional scheduler that attains this maximum [16]. The same holds for minimal reachability probabilities and maximal or minimal expected rewards.
Example 2. Consider the MDP M from Fig. 1(b). We are interested in the maximal probability to reach state s 4 given by Pr max (♦{s 4 }). Since s 4 is not reachable from s 3 we have Pr max s3 (♦{s 4 }) = 0. Intuitively, choosing action β at state s 0 makes reaching s 3 more likely, which should be avoided in order to maximize the probability to reach s 4 . We therefore assume a scheduler σ that always chooses action α at state s 0 . Starting from the initial state s 0 , we then eventually take the transition from s 2 to s 3 or the transition from s 2 to s 4 with probability one. The resulting probability to reach s 4 is given by Pr max (♦{s 4 }) = Pr σ (♦{s 4 }) = 0.3/(0.1 + 0.3) = 0.75.

Probabilistic Model Checking via Interval Iteration
In the following we present approaches to compute reachability probabilities and expected rewards. We consider approximative computations. Exact computations are handled in e.g. [17,18] For the sake of clarity, we focus on reachability probabilities and sketch how the techniques can be lifted to expected rewards.
Reachability Probabilities. We fix an MDP M = (S, Act , P, s I , ρ), a set of goal states G ⊆ S, and a precision parameter ε > 0.
We briefly sketch how to compute such a value r via interval iteration [12,13,19]. The computation for minimal reachability probabilities is analogous.
W.l.o.g. it is assumed that the states in G are absorbing. Using graph algorithms, we compute S 0 = {s ∈ S | Pr max s (♦G) = 0} and partition the state space We say that M is contracting with respect to S ′ ⊆ S if Pr σ s (♦S ′ ) = 1 for all s ∈ S and for all σ ∈ S M . We assume that M is contracting with respect to G ∪ S 0 . Otherwise, we apply a transformation on the so-called end components 1 of M, yielding a contracting MDP M ′ with the same maximal reachability probability as M. Roughly, this transformation replaces each end component of M with a single state whose enabled actions coincide with the actions that previously lead outside of the end component. This step is detailed in [13,19].
We have x * [s] = Pr max s (♦G) for s ∈ S and the unique fixpoint x * of the function f : for s ∈ S ? . Hence, computing Pr max (♦G) reduces to finding the fixpoint of f . A popular technique for this purpose is the value iteration algorithm [1]. Given a starting vector holds for a predefined precision ε > 0. As pointed out in, e.g., [13], there is no guarantee on the preciseness of the result r = f k (x)[s I ], i.e., standard value iteration does not give any evidence on the error |r − Pr max (♦G)|. The intuitive reason is that value iteration only approximates the fixpoint x * from one side, yielding no indication on the distance between the current result and x * . Example 3. Consider the MDP M from Fig. 1(b). We invoked standard value iteration in PRISM [7] and Storm [8] to compute the reachability probability Pr max (♦{s 4 }). Recall from Example 2 that the correct solution is 0.75. With (absolute) precision ε = 10 −6 both model checkers returned 0.7248. Notice that the user can improve the precision by considering, e.g., ε = 10 −8 which yields 0.7497. However, there is no guarantee on the preciseness of a given result.
The interval iteration algorithm [13,12,19] addresses the impreciseness of value iteration. The idea is to approach the fixpoint x * from below and from above. The first step is to find starting vectors we obtain that |r − Pr max (♦G)| < ε, i.e., we get a sound approximation of the maximal reachability probability. Expected Rewards. Whereas [13,19] only consider reachability probabilities, [12] extends interval iteration to compute expected rewards. Let M be an MDP and G be a set of absorbing states such that M is contracting with respect to G.

Problem 2.
Compute an ε-approximation of the maximal expected reachability reward E max ( G), i.e., compute a value r ∈ R with |r − E max ( G)| < ε.
We have for s / ∈ G. As for reachability probabilities, interval iteration can be applied to approximate this fixpoint. The crux lies in finding appropriate starting vectors To this end, [12] describe graph based algorithms that give an upper bound on the expected number of times each individual state s ∈ S \ G is visited. This then yields an approximation of the expected amount of reward collected at the various states.

Sound Value Iteration for MCs
We present an algorithm for computing reachability probabilities and expected rewards as in Problems 1 and 2. The algorithm is an alternative to the interval iteration approach [20,12] but (i) does not require an a priori computation of starting vectors x ℓ , x u ∈ R |S| and (ii) converges faster on many practical benchmarks as shown in Section 5. For the sake of simplicity, we first restrict to computing reachability probabilities on MCs.
In the following, let D = (S, P, s I , ρ) be an MC, G ⊆ S be a set of absorbing goal states and ε > 0 be a precision parameter. We consider the partition S = S 0 ∪ · G ∪ · S ? as in Section 2.2. The following theorem captures the key insight of our algorithm.
Theorem 1. For MC D let G and S ? be as above and k ≥ 0 with Pr s ( ≤k S ? ) < 1 for all s ∈ S ? . We have .
Theorem 1 allows us to approximate Pr(♦G) by computing for increasing k ∈ N -Pr(♦ ≤k G), the probability to reach a state in G within k steps, and -Pr( ≤k S ? ), the probability to stay in S ? during the first k steps. This can be realized via a value-iteration based procedure. The obtained bounds on Pr(♦G) can be tightened arbitrarily since Pr( ≤k S ? ) approaches 0 for increasing k. In the following, we address the correctness of Theorem 1, describe the details of our algorithm, and indicate how the results can be lifted to expected rewards.

Approximating Reachability Probabilities
To approximate the reachability probability Pr(♦G), we consider the step bounded reachability probability Pr(♦ ≤k G) for k ≥ 0 and provide a lower and an upper bound for the 'missing' probability Pr(♦G) − Pr(♦ ≤k G). Note that ♦G is the disjoint union of the paths that reach G within k steps (given by ♦ ≤k G) and the paths that reach G only after k steps (given by S ? U >k G).
A path π ∈ S ? U >k G reaches some state s ∈ S ? after exactly k steps. This yields Consider ℓ, u ∈ [0, 1] with ℓ ≤ Pr s (♦G) ≤ u for all s ∈ S ? , i.e., ℓ and u are lower and upper bounds for the reachability probabilities within S ? . We have s∈S ?
We can argue similar for the lower bound ℓ. With Lemma 1 we get the following. Proposition 1. For MC D with G, S ? , ℓ, u as above and any k ≥ 0 we have We now discuss how the bounds ℓ and u can be obtained from the step bounded probabilities Pr s (♦ ≤k G) and Pr s ( ≤k S ? ) for s ∈ S ? . We focus on the upper bound u. The reasoning for the lower bound ℓ is similar.
Let s max ∈ S ? be a state with maximal reachability probability, that is s max ∈ arg max s∈S ? Pr s (♦G). From Proposition 1 we get We solve the inequality for Pr smax (♦G) (assuming Pr s ( ≤k S ? ) < 1 for all s ∈ S ? ): .
Proposition 2. For MC D let G and S ? be as above and k ≥ 0 such that Pr s ( ≤k S ? ) < 1 for all s ∈ S ? . For everyŝ ∈ S ? we have min s∈S ? .
Theorem 1 is a direct consequence of Propositions 1 and 2.

Extending the Value Iteration Approach
Recall the standard value iteration algorithm for approximating Pr(♦G) as discussed in Section 2. Algorithm 1 depicts our approach. It maintains vectors x k , y k ∈ R |S| which, after k iterations of the loop, store the k-step bounded probabilities Pr s (♦ ≤k G) and Pr s ( ≤k S ? ), respectively. Additionally, the algorithm considers lower bounds ℓ k and upper bounds u k such that the following invariant holds. Input : MC D = (S, P, sI , ρ), absorbing states G ⊆ S, precision ε > 0 Output : Algorithm 1: Sound value iteration for MCs.
Example 5. We apply Algorithm 1 for the MC in Fig. 1(a)

Sound Value Iteration for Expected Rewards
We lift our approach to expected rewards in a straightforward manner. Let G ⊆ S be a set of absorbing goal states of MC D such that Pr(♦G) = 1. Further let S ? = S \ G. For k ≥ 0 we observe that the expected reward E( G) can be split into the expected reward collected within k steps and the expected reward collected only after k steps, i.e., E( G) = E( ≤k G) + s∈S ? Pr(S ? U =k {s}) · E s ( G). Following a similar reasoning as in Section 3.1 we can show the following.
Theorem 3. For MC D let G and S ? be as before and k ≥ 0 such that Pr s ( ≤k S ? ) < 1 for all s ∈ S ? . We have .
Recall the function g : R |S| → R |S| from Section 2.2, given by g(x)[G] = 0 and We modify Algorithm 1 such that it considers function g instead of function f . Then, the returned value r satisfies |r − E( G)| < ε.

Optimizations.
Algorithm 1 can make use of initial bounds ℓ 0 , u 0 ∈ R with ℓ 0 ≤ Pr s (♦G) ≤ u 0 for all s ∈ S ? . Such bounds could be derived, e.g., from domain knowledge or during preprocessing [12]. The algorithm always chooses the largest available lower bound for ℓ k and the lowest available upper bound for u k , respectively. If Algorithm 1 and interval iteration are initialized with the same bounds, Algorithm 1 always requires as most as many iterations compared to interval iteration (cf.

Remark 1).
Gauss-Seidel value iteration [1,12] is an optimization for standard value iteration and interval iteration that potentially leads to faster convergence. When computing f (x)[s] for s ∈ S ? , the idea is to consider already computed results f (x)[s ′ ] from the current iteration. Formally, let ≺ ⊆ S × S be some strict total ordering of the states. Gauss-Seidel value iteration considers instead of function f the function f ≺ : Values f ≺ (x)[s] for s ∈ S are computed in the order defined by ≺. This idea can also be applied to our approach. To this end, we replace f by f ≺ and h by h ≺ , where h ≺ is defined similarly. More details are given in Appendix A.
Topological value iteration [14] employs the graphical structure of the MC D. The idea is to decompose the states S of D into strongly connected components 2 (SCCs) that are analyzed individually. The procedure can improve the runtime of classical value iteration since for a single iteration only the values for the current SCC have to be updated. A topological variant of interval iteration is introduced in [12]. Given these results, sound value iteration can be extended similarly.

Sound Value Iteration for MDPs
We extend sound value iteration to compute reachability probabilities in MDPs. Assume an MDP M = (S, Act, P, s I , ρ) and a set of absorbing goal states G. For simplicity, we focus on maximal reachability probabilities, i.e., we compute Pr max (♦G). Minimal reachability probabilities and expected rewards are analogous. As in Section 2.2 we consider the partition S = S 0 ∪ · G ∪ · S ? such that M is contracting with respect to G ∪ S 0 .

Approximating Maximal Reachability Probabilities
We argue that our results for MCs also hold for MDPs under a given scheduler σ ∈ S M . Let k ≥ 0 such that Pr σ s ( ≤k S ? ) < 1 for all s ∈ S ? . Following the reasoning as in Section 3.1 we get Next, assume an upper bound u ∈ R with Pr max s (♦G) ≤ u for all s ∈ S ? . For a scheduler σ max ∈ S M that attains the maximal reachability probability, i.e., σ max ∈ arg max σ∈S M Pr σ (♦G) it holds that We obtain the following theorem which is the basis of our algorithm.
Theorem 4. For MDP M let G, S ? , and u be as above. Assume σ ∈ S M and k ≥ 0 such that σ ∈ arg max σ ′ ∈S M Pr σ ′ (♦ ≤k G) + Pr σ ′ ( ≤k S ? ) · u and Pr σ s ( ≤k S ? ) < 1 for all s ∈ S ? . We have Similar to the results for MCs it also holds that Pr max (♦G) ≤ max σ∈S Mû σ k witĥ .
However, this upper bound can not trivially be embedded in a value iteration based procedure. Intuitively, in order to compute the upper bound for iteration k, one can not necessarily build on the results for iteration k − 1.

Extending the Value Iteration Approach
The idea of our algorithm is to compute the bounds for Pr max (♦G) as in Theorem 4 for increasing k ≥ 0. Algorithm 2 outlines the procedure. Similar to Algorithm 1 for MCs, vectors x k , y k ∈ R |S| store the step bounded probabilities Pr σ k s (♦ ≤k G) and Pr σ k s ( ≤k S ? ) for any s ∈ S. In addition, schedulers σ k and upper bounds u k ≥ max s∈S ? Pr max s (♦G) are computed in a way that Theorem 4 is applicable. (♦G) ≤ u k , where σ k ∈ arg max σ∈S M Pr σ s (♦ ≤k G) + Pr σ s ( ≤k S ? ) · u k . The lemma holds for k = 0 as x 0 , y 0 , and u 0 are initialized accordingly. For k > 0 we assume that the claim holds after k − 1 iterations, i.e., for x k−1 , y k−1 , u k−1 and scheduler σ k−1 . The results of the kth iteration are obtained as follows.
The function findAction illustrated in Algorithm 3 determines the choices of a scheduler σ k ∈ arg max σ∈S M Pr σ s (♦ ≤k G) + Pr σ s ( ≤k S ? ) · u k−1 for s ∈ S ? . The idea is to consider at state s an action σ k (s) = α ∈ Act(s) that maximizes For the case where no real upper bound is known (i.e., u k−1 = ∞) we implicitly assume a sufficiently large value for u k−1 such that Pr σ s (♦ ≤k G) becomes negligible. Upon leaving state s, σ k mimics σ k−1 , i.e., we set σ k (sαs 1 α 1 . . . s n ) = σ k−1 (s 1 α 1 . . . s n ). After executing Line 15 of Algorithm 2 we have x k [s] = Pr σ k s (♦ ≤k G) and y k [s] = Pr σ k s ( ≤k S ? ). It remains to derive an upper bound u k . To ensure that Lemma 3 holds we require (i) u k ≥ max s∈S ? Pr max s (♦G) and (ii) u k ∈ U k , where Input : MDP M = (S, Act , P, sI , ρ), absorbing states G ⊆ S, precision ε > 0 Output :

Algorithm 2: Sound value iteration for MDPs
Intuitively, the set U k ⊆ R consists of all possible upper bounds u for which σ k is still optimal. U k ⊆ is convex as it can be represented as a conjunction of inequalities with U 0 = R and u ∈ U k if and only if u ∈ U k−1 and for all s ∈ S ? with σ k (s) = α and for all β ∈ Act(s) \ {α} The algorithm maintains the so-called decision value d k which corresponds to the minimum of U k (or −∞ if the minimum does not exist). Algorithm 4 outlines the procedure to obtain the decision value at a given state. Our algorithm ensures that u k is only set to a value in [d k , u k−1 ] ⊆ U k . To show that u k is a valid upper bound, let s max ∈ arg max s∈S ? Pr max s (♦G) and u * = Pr max smax (♦G). From Theorem 4, u k−1 ≥ u * , and u k−1 ∈ U k we get Algorithm 3: Computation of optimal action.
. We distinguish three cases to show that u k = min(u k−1 , max(d k , max s∈S ? Example 7. Reconsider the MDP M from Fig. 2(a) and goal states G = {s 3 , s 4 }.
The maximal reachability probability is attained for a scheduler that always chooses β at state s 0 , which results in Pr max (♦G) = 0.5. We now illustrate how Algorithm 2 approximates this value by sketching the first two iterations. For the first iteration findAction yields action α at s 0 . We obtain: In the second iteration findAction yields again α for s 0 and we get:  Due to the decision value we do not set the upper bound u 2 to 0.29 < Pr max (♦G).
The correctness of the algorithm follows from Theorem 4 and Lemma 3. Termination follows since M is contracting with respect to S 0 ∪ G, implying lim k→∞ Pr σ ( ≤k S ? ) = 0. The optimizations for Algorithm 1 mentioned in Section 3.4 can be applied to Algorithm 2 as well.

Experimental Evaluation
Implementation. We implemented sound value iteration for MCs and MDPs into the model checker Storm [8]. The implementation computes reachability probabilities and expected rewards using explicit data structures such as sparse matrices and vectors. Moreover, Multi-objective model checking is supported, where we straightforwardly extend the value iteration-based approach of [21] to sound value iteration. We also implemented the optimizations given in Section 3. all MCs, MDPs, and CTMCs from the PRISM benchmark suite [22], -several case studies from the PRISM website www.prismmodelchecker.org, -Markov automata accompanying IMCA [23], and multi-objective MDPs considered in [21]. In total, 130 model and property instances were considered. For CTMCs and Markov automata we computed (untimed) reachability probabilities or expected rewards on the underlying MC and the underlying MDP, respectively. In all experiments the precision parameter was given by ε = 10 −6 .
We compare sound value iteration (SVI) with interval iteration (II) as presented in [13,12]. We consider the Gauss-Seidel variant of the approaches and compute initial bounds ℓ 0 and u 0 as in [12]. For a better comparison we consider the implementation of II in Storm. Appendix B gives a comparison with the implementation of II in PRISM. The experiments were run on a single core (2GHz) of an HP BL685C G7 with 192GB of available memory. However, almost all experiments required less than 4GB. We measured model checking times and required iterations. All logfiles and considered benchmarks are available at [24]. Fig. 3(a) depicts the model checking times for SVI (x-axis) and II (y-axis). For better readability, the benchmarks are divided into four plots with different scales. Triangles ( ) and circles (•) indicate MC and MDP benchmarks, respectively. Similarly, Fig. 3(b) shows the required iterations of the approaches. We observe that SVI converged faster and required fewer iterations for almost all MCs and MDPs. SVI performed particularly well on the challenging instances where many iterations are required. Similar observations were made when comparing the topological variants of SVI and II. Both approaches were still competitive if no a priori bounds are given to SVI. More details are given in Appendix B. Fig. 4 indicates the model checking times of SVI and II as well as their topological variants. For reference, we also consider standard (unsound) value iteration (VI). The x-axis depicts the number of instances that have been solved by the corresponding approach within the time limit indicated on the y-axis. Hence, a point (x, y) means that for x instances the model checking time was less or equal than y. We observe that the topological variant of SVI yielded the best run times among all sound approaches and even competes with (unsound) VI.

Conclusion
In this paper we presented a sound variant of the value iteration algorithm which safely approximates reachability probabilities and expected rewards in MCs and MDPs. Experiments on a large set of benchmarks indicate that our approach is a reasonable alternative to the recently proposed interval iteration algorithm.

A Correctness of Sound Gauss-Seidel Value Iteration
We provide additional details on the correctness of the Gauss-Seidel variant of our approach as presented in Section 3.2. Let D = (S, P, s I , ρ) be a DTMC with absorbing set of goal states G and partition S = S 0 ∪ · G ∪ · S ? . Further let ≺ ∈ S × S be an arbitrary strict total order on the states.
Intuitively, the paths in ♦ ≤k κ G reach G within k steps, where only steps from s to s ′ with s ≺ s ′ are counted.
The correctness of the Gauss-Seidel variant of Algorithm 1 follows from the following theorem which is analogue to Theorem 1 and can be shown in a similar way. .

B Additional Experimental Results
We provide additional results of our experiments in order to substantiate our claims from Section 5.   6. Comparison of SVI without computation of a priori bounds (x-axis) and II (y-axis).
In Fig. 5 a comparison of SVI with the implementation of II in PRISM is given. We consider PRISM 4.4 available on its website www.prismmodelchecker.org. Notice that we only consider a subset of our benchmark instances that are supported by both tools. In particular, PRISM currently does not support Markov automata or interval iteration for multi-objective queries.
In Fig. 6 we compare II with a variant of SVI for which no initial bounds u 0 , ℓ 0 were computed.
The topological variants of SVI and II are compared in Fig. 7.