key: cord-0059559-jgl7kuf2
authors: Meyer, Fabian; Hark, Marcel; Giesl, Jürgen
title: Inferring Expected Runtimes of Probabilistic Integer Programs Using Expected Sizes
date: 2021-03-01
journal: Tools and Algorithms for the Construction and Analysis of Systems
DOI: 10.1007/978-3-030-72016-2_14
sha: fd9bec8ddd5928de5a680ce7f849982a2782a7b0
doc_id: 59559
cord_uid: jgl7kuf2

We present a novel modular approach to infer upper bounds on the expected runtimes of probabilistic integer programs automatically. To this end, it computes bounds on the runtimes of program parts and on the sizes of their variables in an alternating way. To evaluate its power, we implemented our approach in a new version of our open-source tool KoAT.

There exist several approaches and tools for automatic complexity analysis of non-probabilistic programs, e.g., [2-6, 8, 9, 18, 20, 21, 27, 28, 30, 34-36, 51, 57, 58] . While most of them rely on basic techniques like ranking functions (see, e.g., [6, 12-14, 17, 53] ), they usually combine these basic techniques in sophisticated ways. For example, in [18] we developed a modular approach for automated complexity analysis of integer programs, based on an alternation between finding symbolic runtime bounds for program parts and using them to infer bounds on the sizes of variables in such parts. So each analysis step is restricted to a small part of the program. The implementation of this approach in KoAT [18] (which is integrated in AProVE [30] ) is one of the leading tools for complexity analysis [31] .

While there exist several adaptions of basic techniques like ranking functions to probabilistic programs (e.g., [1, 11, 15, 16, 22-26, 29, 32, 37, 38, 48, 62] ), most of the sophisticated full approaches for complexity analysis have not been adapted to probabilistic programs yet, and there are only few powerful tools available which analyze the runtimes of probabilistic programs automatically [10, 50, 61, 62] .

We study probabilistic integer programs (Sect. 2) and define suitable notions of non-probabilistic and expected runtime and size bounds (Sect. 3). Then, we adapt our modular approach for runtime and size analysis of [18] to probabilistic programs (Sect. 4 and 5) . So such an adaption is not only possible for basic techniques like ranking functions, but also for full approaches for complexity analysis.

For this adaption, several problems had to be solved. When computing expected runtime or size bounds for new program parts, the main difficulty is to determine when it is sound to use expected bounds on previous program parts and when one has to use non-probabilistic bounds instead. Moreover, the semantics of probabilistic programs is significantly different from classical integer programs. Thus, the proofs of our techniques differ substantially from the ones in [18] , e.g., we have to use concepts from measure theory like ranking supermartingales.

In Sect. 6, we evaluate the implementation of our new approach in the tool KoAT [18, 43] and compare with related work. We refer to [47] for an appendix of our paper containing all proofs, preliminaries from probability and measure theory, and an overview on the benchmark collection used in our evaluation.

For any set M ⊆ R (with R = R ∪ {∞}) and w ∈ M , let M ≥w = {v ∈ M | v ≥ w ∨ v = ∞}. For a set PV of program variables, we first introduce the kind of bounds that our approach computes. Similar to [18] , our bounds represent weakly monotonically increasing functions from PV → R ≥0 . Such bounds have the advantage that they can easily be "composed", i.e., if f and g are both weakly monotonically increasing upper bounds, then so is f • g.

Our notion of probabilistic programs combines classical integer programs (as in, e.g., [18] ) and probabilistic control flow graphs (see, e.g., [1] ). A state s is a variable assignment s : V → Z for the (finite) set V of all variables in the program, where PV ⊆ V, V \ PV is the set of temporary variables, and Σ is the set of all states. For any s ∈ Σ, the state |s| is defined by |s| (x) = |s(x)| for all x ∈ V. The set C of constraints is the smallest set containing e 1 ≤ e 2 for all polynomials e 1 , e 2 ∈ Z[V] and c 1 ∧ c 2 for all c 1 , c 2 ∈ C. In addition to "≤", in examples we also use relations like ">", which can be simulated by constraints (e.g., e 1 > e 2 is equivalent to e 2 + 1 ≤ e 1 when regarding integers). We also allow the application of states to arithmetic expressions e and constraints c. Then the number s(e) resp. s(c) ∈ {t, f } results from evaluating the expression resp. the constraint when substituting every variable x by s(x). So for bounds b ∈ B, we have |s| (b) ∈ R ≥0 .

In the transitions of a program, a program variable x ∈ PV can also be updated by adding a value according to a bounded distribution function d : Σ → Dist(Z). Here, for any state s, d(s) is the probability distribution of the values that are added to x. As usual, a probability distribution on Z is a mapping pr : Z → R with pr(v) ∈ [0, 1] for all v ∈ Z and v∈Z pr(v) = 1. Let Dist(Z) be the set of distributions pr whose expected value E(pr) = v∈Z v · pr(v) is well defined and finite, i.e., E abs (pr)

Let D denote the set of all bounded distribution functions (our implementation supports Bernoulli, uniform, geometric, hypergeometric, and binomial distributions, see [43] for details). ∪ D of t, mapping every program variable to an update polynomial or a bounded distribution function. All t ∈ g must have the same start location and the same guard τ . Thus, we call them the start location and guard of g, and denote them by g and τ g . Moreover, the probabilities p of the transitions in g must add up to 1. 4 . an initial location 0 ∈ L, where no transition has target location 0 PIPs allow for both probabilistic and non-deterministic branching and sampling. Probabilistic branching is modeled by selecting a transition out of a non-singleton general transition. Non-deterministic branching is represented by several general transitions with the same start location and non-exclusive guards. Probabilistic sampling is realized by update functions that map a program variable to a bounded distribution function. Non-deterministic sampling is modeled by updating a program variable with an expression containing temporary variables from V \ PV, whose values are non-deterministic (but can be restricted in the guard). The set of initial general transitions GT 0 ⊆ GT consists of all general transitions with start location 0 . Fig. 1 with initial location 0 and the program variables PV = {x, y}. Here, let p = 1 and τ = t if not stated explicitly. There are four general transitions: g 0 = {t 0 }, g 1 = {t 1 , t 2 }, g 2 = {t 3 }, and g 3 = {t 4 }, where g 1 and g 2 represent a non-deterministic branching. When choosing the general transition g 1 , the transitions t 1 and t 2 encode a probabilistic branching. If we modified the update η and the guard τ of t 0 to η(x) = u ∈ V \PV and τ = (u > 0), then x would be updated to a non-deterministically chosen positive value. In contrast, if η(x) = GEO( 1 2 ), then t 0 would update x by adding a value sampled from the geometric distribution with parameter 1 2 . In the following, we regard a fixed PIP P as in Def. 2. A configuration is a tuple ( , t, s), with the current location ∈ L, the current state s ∈ Σ, and the transition t that was evaluated last and led to the current configuration. Let T = g∈GT g.

×Σ is the set of all configurations, with a special location ⊥ indicating the termination of a run, and special transitions t in (used in the first configuration of a run) and t ⊥ (for the configurations of the run after termination). The (virtual) general transition g ⊥ = {t ⊥ } only contains t ⊥ .

A run of a PIP is an infinite sequence ϑ = c 0 c 1 · · · ∈ Conf ω . Let Runs = Conf ω and let FPath = Conf * be the set of all finite paths of configurations.

In our setting, deterministic Markovian schedulers suffice to resolve all nondeterminism (see, e.g., [54, Prop. 6.2.1] ). For c = ( , t, s) ∈ Conf, such a scheduler S yields a pair S(c) = (g, s ) where g is the next general transition to be taken (with = g ) and s chooses values for the temporary variables where s (τ g ) = t and s(x) = s (x) for all x ∈ PV. If GT contains no such g, we get S(c) = (g ⊥ , s).

For each scheduler S and initial state s 0 , we first define a probability mass function pr S,s0 . For all c ∈ Conf, pr S,s0 (c) is the probability that a run starts in c. Thus, pr S,s0 (c) = 1 if c = ( 0 , t in , s 0 ) and pr S,s0 (c) = 0, otherwise. Moreover, for all c , c ∈ Conf, pr S,s0 (c → c) is the probability that the configuration c is followed by the configuration c (see [47] for the formal definition of pr S,s0 ).

For any f = c 0 · · · c n ∈ FPath, let pr S,s0 (f ) = pr S,s0 (c 0 ) · pr S,s0 (c 0 → c 1 ) · . . . · pr S,s0 (c n−1 → c n ). We say that f is admissible for S and s 0 if pr S,s0 (f ) > 0. A run ϑ is admissible if all its finite prefixes are admissible. A configuration c ∈ Conf is admissible if there is some admissible finite path ending in c.

The semantics of PIPs can now be defined by giving a corresponding probability space, which is obtained by a standard cylinder construction (see, e.g., [7, 60] ). Let P S,s0 denote the corresponding probability measure which lifts pr S,s0 to cylinder sets: For any f ∈ FPath, we have pr S,s0 (f ) = P S,s0 (Pre f ) for the set Pre f of all runs with prefix f . So P S,s0 (Θ) is the probability that a run from Θ ⊆ Runs is obtained when using the scheduler S and starting in s 0 .

We denote the associated expected value operator by E S,s0 . So for any random variable X : Runs → N = N ∪ {∞}, we have E S,s0 (X) = n∈N n · P S,s0 (X = n). For details on the preliminaries from probability theory we refer to [47] .

In Sect. 3.1, we first recapitulate the concepts of (non-probabilistic) runtime and size bounds from [18] . Then we introduce expected runtime and size bounds in Sect. 3.2 and connect them to their non-probabilistic counterparts.

Again, let P denote the PIP which we want to analyze. Def. 4 recapitulates the notions of runtime and size bounds from [18] in our setting. Recall that bounds from B do not contain temporary variables, i.e., we always try to infer bounds in terms of the initial values of the program variables. Let sup ∅ = 0, as all occurring sets are subsets of R ≥0 , whose minimal element is 0. [18] ). RB : T → B is a runtime bound and SB : T × V → B is a size bound if for all transitions t ∈ T , all variables x ∈ V, all schedulers S, and all states s 0 ∈ Σ, we have

So RB(t) is a bound on the number of executions of t and SB(t, x) overapproximates the greatest absolute value that x ∈ V takes after the application of the transition t in any admissible finite path. Note that Def. 4 does not apply to t in and t ⊥ , since they are not contained in T .

We call a tuple (RB, SB) a (non-probabilistic) bound pair. We will use such non-probabilistic bound pairs for an initialization of expected bounds (Thm. 10) and to compute improved expected runtime and size bounds in Sect. 4 and 5.

Example 5 (Bound Pair). The technique of [18] computes the following bound pair for the PIP of Fig. 1 (by ignoring the probabilities of the transitions).

Clearly, t 0 and t 3 can only be evaluated once. Since t 1 decrements x and no transition increments it, t 1 's runtime is bounded by |s 0 | (x). However, t 2 can be executed arbitrarily often if s 0 (x) > 0. Thus, the runtimes of t 2 and t 4 are unbounded (i.e., P is not terminating when regarding it as a non-probabilistic program). SB(t, x) is finite for all transitions t, since x is never increased. In contrast, the value of y can be arbitrarily large after all transitions but t 0 .

We now define the expected runtime and size complexity of a PIP P.

Definition 6 (Expected Runtime Complexity, PAST [15] ). For g ∈ GT , its runtime is the random variable R(g) where R : GT → Runs → N with

For a scheduler S and s 0 ∈ Σ, the expected runtime complexity of g ∈ GT is E S,s0 (R(g)) and the expected runtime complexity of P is g∈GT E S,s0 (R(g)).

If P's expected runtime complexity is finite for every scheduler S and every initial state s 0 , then P is called positively almost surely terminating (PAST).

So R(g)(ϑ) is the number of executions of a transition from g in the run ϑ.

While non-probabilistic size bounds refer to pairs (t, x) of transitions t ∈ T and variables x ∈ V (so-called result variables in [18] ), we now introduce expected size bounds for general result variables (g, , x), which consist of a general transition g, one of its target locations , and a program variable x ∈ PV. So x must not be a temporary variable (which represents non-probabilistic non-determinism), since general result variables are used for expected size bounds.

For a scheduler S and s 0 , the expected size complexity of α ∈GRV is E S,s0 (S(α)).

So for any run ϑ, S(g, , x)(ϑ) is the greatest absolute value of x in location , whenever was entered with a transition from g. We now define bounds for the expected runtime and size complexity which hold independent of the scheduler.

• RB E : GT → B is an expected runtime bound if for all g ∈ GT , all schedulers S, and all s 0 ∈ Σ, we have |s 0 | (RB E (g)) ≥ E S,s0 (R(g)).

Sect. 4 and 5 will derive the following expected bounds for the PIP from Fig. 1 .

While the runtimes of t 2 and t 4 were unbounded in the non-probabilistic case (Ex. 5), we obtain finite bounds on the expected runtimes of

For example, we can expect x to be non-positive after at most |s 0 | (2·x) iterations of g 1 . Based on the above expected runtime bounds, the expected runtime complexity of the PIP is at most

where n is the maximal absolute value of the program variables at the start of the program.

The following theorem shows that non-probabilistic bounds can be lifted to expected bounds, since they do not only bound the expected value of R(g) resp. S(α), but the whole distribution. As mentioned, all proofs can be found in [47] .

Here, we over-approximate the maximum of SB(t, x) for t = ( , , , , ) ∈ g by their sum. For asymptotic bounds, this does not affect precision, since max(f, g) and f + g have the same asymptotic growth for any non-negative functions f, g.

Example 11 (Lifting of Bounds). When lifting the bound pair of Ex. 5 to expected bounds according to Thm. 10, one would obtain RB E (g 0 ) = RB E (g 2 ) = 1 and 1 , y) = y, and SB E (g, , y) = ∞ whenever g = g 0 . Thus, with these lifted bounds one cannot show that P's expected runtime complexity is finite, i.e., they are substantially less precise than the finite expected bounds from Ex. 9. Our approach will compute such finite expected bounds by repeatedly improving the lifted bounds of Thm. 10.

We first present a new variant of probabilistic linear ranking functions in Sect. 4.1. Based on this, in Sect. 4.2 we introduce our modular technique to infer expected runtime bounds by using expected size bounds.

For probabilistic programs, several techniques based on ranking supermartingales have been developed. In this section, we define a class of probabilistic ranking functions that will be suitable for our modular analysis.

We restrict ourselves to ranking functions r : L → R[PV] lin that map every location to a linear polynomial (i.e., of at most degree 1) without temporary variables. The linearity restriction is common to ease the automated inference of ranking functions. Moreover, this restriction will be needed for the soundness of our technique. Nevertheless, our approach of course also infers non-linear expected runtimes (by combining the linear bounds obtained for different program parts).

Let exp r,g,s denote the expected value of r after an execution of g ∈ GT in state s ∈ Σ. Here, s η (x) is the expected value of x ∈ PV after performing the update η in state s. So if η(x) ∈ D, then x's expected value after the update results from adding the expected value of the probability distribution η(x)(s):

lin is a probabilistic linear ranking function (PLRF) for GT > and GT ni if for all g ∈ GT ni \ GT > and c ∈ Conf there is a g,c ∈ {<, ≥} such that for all finite paths · · · c c that are admissible for some S and s 0 ∈ Σ, and where c = ( , t, s) (i.e., where t is the transition that is used in the step from c to c), we have: Boundedness (a): If t ∈ g for a g ∈ GT ni \ GT > , then s(r( )) g,c 0. Boundedness (b): If t ∈ g for a g ∈ GT > , then s(r( )) ≥ 0. Non-Increase: If = g for a g ∈ GT ni and s(τ g ) = t, then s(r( )) ≥ exp r,g,s . Decrease: If = g for a g ∈ GT > and s(τ g ) = t, then s(r( )) − 1 ≥ exp r,g,s .

So if one is restricted to the sub-program with the non-increasing transitions GT ni , then r( ) is an upper bound on the expected number of applications of transitions from GT > when starting in . Hence, a PLRF for GT > = GT ni = GT would imply that the program is PAST (see, e.g., [1, 16, 24, 25] ). However, our PLRFs differ from the standard notion of probabilistic ranking functions by considering arbitrary subsets GT ni ⊆ GT . This is needed for the modularity of our approach which allows us to analyze program parts separately (e.g., GT \GT ni is ignored when inferring a PLRF). Thus, our "Boundedness" conditions differ slightly from the corresponding conditions in other definitions. Condition (b) requires that g ∈ GT > never leads to a configuration where r is negative. Condition (a) states that in an admissible path where g = {t 1 , t 2 , . . .} ∈ GT ni \ GT > is used for continuing in configuration c , if executing t 1 in c makes r negative, then executing t 2 must make r negative as well. Thus, such a g can never come before a general transition from GT > in an admissible path and hence, g can be ignored when inferring upper bounds on the runtime. This increases the power of our approach and it allows us to consider only non-negative random variables in our correctness proofs.

We use SMT solvers to generate PLRFs automatically. Then for "Boundedness", we regard all s ∈ Σ with s (τ g ) = t and require "Boundedness" for any state s that is reachable from s . Fig. 1 and the sets GT > = GT ni = {g 1 } and GT > = GT ni = {g 3 }, which correspond to its two loops.

The function r with r( 1 ) = 2 · x and r( 0 ) = r( 2 ) = 0 is a PLRF for GT > = GT ni : For every admissible configuration ( , t, s) with t ∈ g 1 we have = 1 and s(r( 1 )) = 2 · s(x) ≥ 0, since x was positive before (due to g 1 's guard) and it was either decreased by 1 or not changed by the update of t 1 resp. t 2 . Hence r is bounded. Moreover, for s 1 (x) = s(x − 1) = s(x) − 1 and s 2 (x) = s(x) we have: ( 1 )) − 1 So r is decreasing on g 1 and as GT > = GT ni , also the non-increase property holds. Similarly, r with r ( 2 ) = y and r ( 0 ) = r ( 1 ) = 0 is a PLRF for GT > = GT ni .

In our implementation, GT > is always a singleton and we let GT ni ⊆ GT be a cycle in the call graph where we find a PLRF for GT > ⊆ GT ni . The next subsection shows how we can then obtain an expected runtime bound for the overall program by searching for suitable ranking functions repeatedly.

Our approach to infer expected runtime bounds is based on an underlying (nonprobabilistic) bound pair (RB, SB) which is computed by existing techniques (in our implementation, we use [18] ). To do so, we abstract the PIP to a standard integer transition system by ignoring the probabilities of transitions and replacing probabilistic with non-deterministic sampling (e.g., the update η(x) = GEO( 1 2 ) would be replaced by η(x) = x + u with u ∈ V \ PV, where u > 0 is added to the guard). Of course, we usually have RB(t) = ∞ for some transitions t.

We start with the expected bound pair (RB E , SB E ) that is obtained by lifting (RB, SB) as in Thm. 10. Afterwards, the expected runtime bound RB E is improved repeatedly by applying the following Thm. 16 (and similarly, SB E is improved repeatedly by applying Thm. 23 and 25 from Sect. 5). Our approach alternates the improvement of RB E and SB E , and it uses expected size bounds on "previous" transitions to improve expected runtime bounds, and vice versa.

To improve RB E , we generate a PLRF r for a part of the program. To obtain a bound for the full program from r, one has to determine which transitions can enter the program part and from which locations it can be entered. Fig. 1 and GT ni = {g 1 }, we have EL GTni = { 1 } and ET GTni ( 1 ) = {g 0 }. So the loop formed by g 1 is entered at location 1 and the general transition g 0 has to be executed before. Similarly, for GT ni = {g 3 

Recall that if r is a PLRF for GT > ⊆ GT ni , then in a program that is restricted to GT ni , r( ) is an upper bound on the expected number of executions of transitions from GT > when starting in . Since r( ) may contain negative coefficients, it is not weakly monotonically increasing in general. To turn expressions e ∈ R[PV] into bounds from B, let the over-approximation · replace all coefficients by their absolute value. So for example, x − y = x + (−1) · y = x + y. Clearly, we have |s| ( e ) ≥ |s| (e) for all s ∈ Σ. Moreover, if e ∈ R[PV] then e ∈ B.

To turn r( ) into a bound for the full program, one has to take into account how often the sub-program with the transitions GT ni is reached via an entry transition h ∈ ET GTni ( ) for some ∈ EL GTni . This can be over-approximated by t=( , , , , )∈h RB(t), which is an upper bound on the number of times that transitions in h to the entry location of GT ni are applied in a full program run.

The bound r( ) is expressed in terms of the program variables at the entry location of GT ni . To obtain a bound in terms of the variables at the start of the program, one has to take into account which value a program variable x may have when the sub-program GT ni is reached. For every entry transition h ∈ ET GTni ( ), this value can be over-approximated by SB E (h, , x) . Thus, we have to instantiate each variable x in r( ) by SB E (h, , x) . Let SB E (h, , ·) : PV → B be the mapping with SB E (h, , ·)(x) = SB E (h, , x). Hence, SB E (h, , ·)( r( ) ) overapproximates the expected number of applications of GT > if GT ni is entered in location , where this bound is expressed in terms of the input variables of the program. Here, weak monotonic increase of r( ) ensures that instantiating its variables by an over-approximation of their size yields an over-approximation of the runtime. Fig. 1 

RB E (g) = ⎧ ⎪ ⎨ ⎪ ⎩ ∈ELGT ni h∈ETGT ni ( ) ( t=( , , , , )∈h RB(t)) · (SB E (h, , ·) ( r( ) )) , if g ∈ GT > RB E (g), if g ∈ GT >

1 For a set of sets like GTni, GTni denotes their union, i.e., GTni = g∈GT ni g.

To improve RB E (g 3 ), we use the PLRF r for GT > = GT ni = {g 3 } from Ex. 13 . As EL GT ni = { 2 } and ET GT ni ( 2 ) = {g 2 } by Ex. 15 , where g 2 = {t 3 } and RB(t 3 ) = 1 (Ex. 5), with the bound SB E (g 2 , 2 , y) = 6 · x 2 + 2 · y from Ex. 9, Thm. 16 yields

So based on the expected size bounds of Ex. 9, we have shown how to compute the expected runtime bounds of Ex. 9 automatically.

Similar to [18] , our approach relies on combining bounds that one has computed earlier in order to derive new bounds. Here, bounds may be combined linearly, bounds may be multiplied, and bounds may even be substituted into other bounds. But in contrast to [18] , sometimes one may combine expected bounds that were computed earlier and sometimes it is only sound to combine non-probabilistic bounds: If a new bound is computed by linear combinations of earlier bounds, then it is sound to use the "expected versions" of these earlier bounds. However, if two bounds are multiplied, then it is in general not sound to use their "expected versions". Thus, it would be unsound to use the expected runtime bounds RB E (h) instead of the non-probabilistic bounds t=( , , , , )∈h RB(t) on the entry transitions in Thm. 16 (a counterexample is given in [47] ). 2 In general, if bounds b 1 , . . . , b n are substituted into another bound b, then it is sound to use "expected versions" of the bounds b 1 , . . . , b n if b is concave, see, e.g., [10, 11, 40] . Since bounds from B do not contain negative coefficients, we obtain that a finite 3 bound b ∈ B is concave iff it is a linear polynomial (see [47] ).

Thus, in Thm. 16 we may substitute expected size bounds SB E (h, , x) into r( ) , since we restricted ourselves to linear ranking functions r and hence, r( ) is also linear. Note that in contrast to [11] , where a notion of concavity was used to analyze probabilistic term rewriting, a multilinear expression like x · y is not concave when regarding both arguments simultaneously. Hence, it is unsound to use such ranking functions in Thm. 16 . See [47] for a counterexample to show why substituting expected bounds into a non-linear bound is incorrect in general.

We first compute local bounds for one application of a transition (Sect. 5.1). To turn them into global bounds, we encode the data flow of a PIP in a graph. Sect. 5.2 then presents our technique to compute expected size bounds.

We first compute a bound on the expected change of a variable during an update. More precisely, for every general result variable (g, , x) we define a bound CB E (g, , x) on the change of the variable x that we can expect in one execution of the general transition g when reaching location . So we consider all t = ( , p, , η, ) ∈ g and the expected difference between the current value of x and its update η(x). However, for η(x) ∈ Z[V], η(x) − x is not necessarily from B because it may contain negative coefficients. Thus, we use the overapproximation η(x) − x (where we always simplify expressions before applying · , e.g., x − x = 0 = 0). Moreover, η(x) − x may contain temporary variables. Let tv t : V → B instantiate all temporary variables by the largest possible value they can have after evaluating the transition t. Hence, we then use tv t ( η(x) − x ) instead. For tv t , we have to use the underlying non-probabilistic size bound SB for the program (since the scheduler determines the values of temporary variables by non-deterministic (non-probabilistic) choice). If x is updated according to a bounded distribution function d ∈ D, then as in Sect. 2, let E(d) ∈ B denote a finite bound on d, i.e., E abs (d(s)) ≤ |s| (E(d)) for all s ∈ Σ.

Then

. For the PIP of Fig. 1 , we have CB E (g 0 , , ) = CB E (g 2 , , ) = CB E (g 3 , 2 , x) = 0, since the respective updates are identities. Moreover,

In a similar way, we obtain CB E (g 1 , 1 , y) = x and CB E (g 3 , 2 , y) = 1. To obtain global bounds from the local bounds CB E (g, , x), we construct a general result variable graph which encodes the data flow between variables. Let pre(g) = ET ∅ ( g ) be the the set of pre-transitions of g which lead into g's start location g . Moreover, for α = (g, , x) ∈ GRV let its active variables actV(α) consist of all variables occurring in the bound x + CB E (α) for α's expected size. Graph) . The general result variable graph has the set of nodes GRV and the set of edges GRVE, where GRVE = { ((g , , x ), (g, , x) ) | g ∈ pre(g) ∧ = g ∧ x ∈ actV(g, , x) }.

The general result variable graph for the PIP of Fig. 1 Similarly, actV(g 1 , 1 , y) = {x, y}, as x and y are contained in y+CB E (g 1 , 1 , y) = y+x. For all other α ∈ GRV, we have actV( , , x) = {x} and actV( , , y) = {y}. As pre(g 1 ) = {g 0 , g 1 }, the graph captures the dependence of (g 1 , 1 , x) on (g 0 , 1 , x) and (g 1 , 1 , x) , and of (g 1 , 1 , y) on (g 0 , 1 , x), (g 0 , 1 , y), (g 1 , 1 , x), and (g 1 , 1 , y) . The other edges are obtained in a similar way.

We now compute global expected size bounds for the general result variables by considering the SCCs of the general result variable graph separately. As usual, an SCC is a maximal subgraph with a path from each node to every other node. An SCC is trivial if it consists of a single node without an edge to itself. We first handle trivial SCCs in Sect. 5.2.1 and consider non-trivial SCCs in Sect. 5.2.2.

By Thm. 20, x + CB E (g, , x) is a local bound on the expected value of x after applying g once in order to enter . However, this bound is formulated in terms of the values of the variables immediately before applying g. We now want to compute global bounds in terms of the initial values of the variables at the start of the program.

If g is initial (i.e., g ∈ GT 0 since g starts in the initial location 0 ), then x + CB E (g, , x) is already a global bound, as the values of the variables before the application of g are the initial values of the variables at the program start.

Otherwise, the variables y occurring in the local bound x + CB E (g, , x) have to be replaced by the values that they can take in a full program run before applying the transition g. Thus, we have to consider all transitions h ∈ pre(g) and instantiate every variable y by the maximum of the values that y can have after applying h. Here, we again over-approximate the maximum by the sum.

If CB E (g, , x) is concave (i.e., a linear polynomial), then we can instantiate its variables by expected size bounds SB E (h, g , y). However, this is unsound if CB E (g, , x) is not linear, i.e., not concave (see [47] for a counterexample). So in this case, we have to use non-probabilistic bounds SB(t, y) instead.

As in Sect. 4.2, we use an underlying non-probabilistic bound pair (RB, SB) and start with the expected pair (RB E , SB E ) obtained by lifting (RB, SB) according to Thm. 10. While Thm. 16 improves RB E , we now improve SB E . Here, the SCCs of the general result variable graph should be treated in topological order, since then one may first improve SB E for result variables corresponding to pre(g), and use that when improving SB E for result variables of the form (g, , ).

Let SB E be an expected size bound, SB a (non-probabilistic) size bound, and let α = (g, , x) form a trivial SCC of the general result variable graph. Let size α E and size α be mappings from PV → B with size α E (y) = h∈pre(g) SB E (h, g , y) and size α (y) = h∈pre(g), t=( , , , , g )∈h SB(t, y). Then SB E : GRV → B is an expected size bound,

where SB E (β) = SB E (β) for β = α and = (g 0 , 1 , x), α y = (g 0 , 1 , y) , β x = (g 2 , 2 , x), and β y = (g 2 , 2 , y). For all these general result variables, the expected local change bound CB E is 0 (see Ex. 19) . Thus, it is linear. Since g 0 ∈ GT 0 , Thm. 23 yields

By treating SCCs in topological order, when handling β x , β y , we can assume that we already have SB E (α x ) = x, SB E (α y ) = y and SB E (g 1 , 1 , x) = 2 · x, SB E (g 1 , 1 , y) = 6 · x 2 + y (see Ex. 9) for the result variables corresponding to pre(g 2 ) = {g 0 , g 1 }. We will explain in Sect. 5 

Now we handle non-trivial SCCs C of the general result variable graph. An upper bound for the expected size of a variable x when entering C is obtained from SB E (β) for all general result variables β = ( , , x) which have an edge to C.

To turn CB E (g, , x) into a global bound, as in Thm. 23 its variables y have to be instantiated by the values size (g, ,x) (y) that they can take in a full program run before applying a transition from g. Thus, size (g, ,x) (CB E (g, , x)) is a global bound on the expected change resulting from one application of g. To obtain an upper bound for the whole SCC C, we add up these global bounds for all (g, , x) ∈ C and take into account how often the general transitions in the SCC are expected to be executed, i.e., we multiply with their expected runtime bound RB E (g). So while in Thm. 16 we improve RB E using expected size bounds for previous transitions, we now improve SB E (C) using expected runtime bounds for the transitions in C and expected size bounds for previous transitions.

Then SB E is an expected size bound:

Here we really have to use the non-probabilistic size bound size α instead of size α E , even if CB E (α ) is linear, i.e., concave. Otherwise we would multiply the expected values of two random variables which are not independent.

Ex. 22 contains 4 non-trivial SCCs formed by α x = (g 1 , 1 , x), α y = (g 1 , 1 , y), β x = (g 3 , 2 , x), and β y = (g 3 , 2 , y). By the results on SB E , RB E , CB E , and SB from Ex. 24, 17, 19, and 5, Thm. 25 yields the expected size bound in Ex. 9 :

Related Work Our approach adapts techniques from [18] to probabilistic programs. As explained in Sect. 1, this adaption is not at all trivial (see our proofs in [47] ). There has been a lot of work on proving PAST and inferring bounds on expected runtimes using supermartingales, e.g., [1, 11, 15, 16, 22-25, 29, 32, 48, 62] . While these techniques infer one (lexicographic) ranking supermartingale to analyze the complete program, our approach deals with information flow between different program parts and analyzes them separately.

There is also work on modular analysis of almost sure termination (AST) [1, 25, 26, 37, 38, 48] , i.e., termination with probability 1. This differs from our results, since AST is compositional, in contrast to PAST (see, e.g., [41, 42] ).

A fundamentally different approach to ranking supermartingales (i.e., to forward-reasoning) is backward-reasoning by so-called expectation transformers, see, e.g., [10, 41, 42, 44-46, 50, 52, 61] . In this orthogonal reasoning, [10, 41, 42, 52] consider the connection of the expected runtime and size. While expectation transformers apply backward-instead of forward-reasoning, their correctness can also be justified using supermartingales. More precisely, Park induction for upper bounds on the expected runtime via expectation transformers essentially ensures that a certain stochastic process is a supermartingale (see [33] for details).

To the best of our knowledge, the only available tools for the inference of upper bounds on the expected runtimes of probabilistic programs are [10, 50, 61, 62] . The tool of [61] deals with data types and higher order functions in probabilistic ML programs and does not support programs whose complexity depends on (possibly negative) integers (see [55] ). Furthermore, the tool of [48] focuses on proving or refuting (P)AST of probabilistic programs for so-called Prob-solvable loops, which do not allow for nested or sequential loops or non-determinism. So both [61] and [48] are orthogonal to our work. We discuss [10, 50, 62] below. Implementation We implemented our analysis in a new version of our tool KoAT [18] .

KoAT is an open-source tool written in OCaml, which can also be downloaded as a Docker image and accessed via a web interface [43] .

Given a PIP, the analysis proceeds as in Alg. 1. The preprocessing in Line 1 adds invariants to guards (using APRON [39] to generate (non-probabilistic) invariants), unfolds transitions [19] , and removes unreachable locations, transitions with probability 0, and transitions with unsatisfiable guards (using Z3 [49] ). Input: PIP (PV, L, GT , 0) 1 preprocess the PIP 2 (RB, SB) ← perform non-probabilistic analysis using [18] 3 (RB E , SB E ) ← lift (RB, SB) to an expected bound pair with Thm. 10 4 repeat 5 for all SCCs C of the general result variable graph in topological order do

for all general transitions g ∈ GT do 10 RB E ← improve RB E for GT> = {g} by Thm. 16 11 RB E (g) ← min{RB E (g), RB E (g)} 12 until no bound is improved anymore Output: g∈GT RB E (g) Algorithm 1: Overall approach to infer bounds on expected runtimes

We start by a non-probabilistic analysis and lift the resulting bounds to an initial expected bound pair (Lines 2 and 3). Afterwards, we first try to improve the expected size bounds using Thm. 23 and 25, and then we attempt to improve the expected runtime bounds using Thm. 16 (if we find a PLRF using Z3). To determine the "minimum" of the previous and the new bound, we use a heuristic which compares polynomial bounds by their degree. While we over-approximated the maximum of expressions by their sum to ease readability in this paper, KoAT also uses bounds containing "min" and "max" to increase precision.

This alternating modular computation of expected size and runtime bounds is repeated so that one can benefit from improved expected runtime bounds when computing expected size bounds and vice versa. We abort this improvement of expected bounds in Alg. 1 if they are all finite (or when reaching a timeout).

To assess the power of our approach, we performed an experimental evaluation of our implementation in KoAT. We did not compare with the tool of [62] , since [62] expects the program to be annotated with already computed invariants. But for many of the examples in our experiments, the invariant generation tool [56] used by [62] did not find invariants strong enough to enable a meaningful analysis (and we could not apply APRON [39] due to the different semantics of invariants).

Instead, we compare KoAT with the tools Absynth [50] and eco-imp [10] which are both based on a conceptionally different backward-reasoning approach. We ran the tools on all 39 examples from Absynth's evaluation in [50] (except recursive, which contains non-tail-recursion and thus cannot be encoded as a PIP), and on the 8 additional examples from the artifact of [50] . Moreover, our collection has 29 additional benchmarks: 14 examples that illustrate different aspects of PIPs, 5 PIPs based on examples from [50] where we removed assumptions, and 10 PIPs based on benchmarks from the TPDB [59] where some transitions were enriched with probabilistic behavior. The TPDB is a collection of typical programs used in the annual Termination and Complexity Competition [31] . We ran the experiments on an iMac with an Intel i5-2500S CPU and 12 GB of RAM under macOS Sierra for Absynth and NixOS 20.03 for KoAT and eco-imp. A timeout of 5 minutes per Fig. 3 : Results on our new benchmarks example was applied for all tools. The average runtime of successful runs was 4.26 s for KoAT, 3.53 s for Absynth, and just 0.93 s for eco-imp. Fig. 2 and 3 show the generated asymptotic bounds, where n is the maximal absolute value of the program variables at the program start. Here, "∞" indicates that no finite time bound could be computed and "TO" means "timeout". The detailed asymptotic results of all tools on all examples can be found in [43, 47] .

Absynth and eco-imp slightly outperform KoAT on the examples from Absynth's collection, while KoAT is considerably stronger than both tools on the additional benchmarks. In particular, Absynth and eco-imp outperform our approach on examples with nested probabilistic loops. While our modular approach can analyze inner loops separately when searching for probabilistic ranking functions, Thm. 16 then requires non-probabilistic time bounds for all transitions entering the inner loop. But these bounds may be infinite if the outer loop has probabilistic behavior itself. Moreover, in contrast to our work and [10] , the approach of [50] does not require weakly monotonic bounds.

On the other hand, KoAT is superior to Absynth and eco-imp on large examples with many loops, where only a few transitions have probabilistic behavior (this might correspond to the typical application of randomization in practical programming). Here, we benefit from the modularity of our approach which treats loops independently and combines their bounds afterwards. Absynth and eco-imp also fail for our leading example of Fig. 1 , while KoAT infers a quadratic bound. Hence, the tools have particular strengths on orthogonal kinds of examples.

KoAT's source code is available at https://github.com/aprove-developers/ KoAT2-Releases/tree/probabilistic. To obtain a KoAT artifact, see https:// aprove-developers.github.io/ExpectedUpperBounds/ for a static binary and Docker image. This web site also provides all examples from our evaluation, detailed outputs of our experiments, and a web interface to run KoAT directly online. Conclusion We presented a new modular approach to infer upper bounds on the expected runtimes of probabilistic integer programs. To this end, non-probabilistic and expected runtime and size bounds on parts of the program are computed in an alternating fashion and then combined to an overall expected runtime bound. In the evaluation, our tool KoAT succeeded on 91% of all examples, while the main other related tools (Absynth and eco-imp) only inferred finite bounds for 68% resp. 77% of the examples. In future work, it would be interesting to consider a modular combination of these tools (resp. of their underlying approaches).

Lexicographic ranking supermartingales: An efficient approach to termination of probabilistic programs

Closed-form upper bounds in static cost analysis

Cost analysis of object-oriented bytecode programs

On the inference of resource usage upper and lower bounds

Resource analysis driven by (conditional) termination proofs. Theory Pract

Multi-dimensional rankings, program termination, and complexity bounds of flowchart programs

Probability and Measure Theory

A combination framework for complexity

TcT: Tyrolean Complexity Tool

A modular cost analysis for probabilistic programs

On probabilistic term rewriting

Ranking functions for linear-constraint loops

On multiphase-linear ranking functions

Multiphase-linear ranking functions and their relation to recurrent sets

Proving positive almost-sure termination

Proving positive almost sure termination under strategies

Linear ranking with reachability

Analyzing runtime and size complexity of integer programs

A transformation system for developing recursive programs

Compositional certified resource bounds

Automated resource analysis with Coq proof objects

Probabilistic program analysis with martingales

Stochastic invariants for probabilistic termination

Algorithmic analysis of qualitative and quantitative termination problems for affine probabilistic programs

Termination analysis of probabilistic programs with martingales

Probabilistic termination: Soundness, completeness, and compositionality

Resource analysis of complex programs with cost equations

Upper and lower amortized cost bounds of programs expressed as cost relations

Termination of nondeterministic probabilistic programs

Analyzing program termination and complexity automatically with AProVE

The termination and complexity competition

Computing expected runtimes for constant probability programs

Aiming low is harder: Induction for lower bounds in probabilistic program verification

Multivariate amortized resource analysis

Type-based amortized resource analysis with integers and arrays

Towards automatic resource bound analysis for OCaml

New approaches for almost-sure termination of probabilistic programs

Modular verification for almost-sure termination of probabilistic programs

Apron: A library of numerical abstract domains for static analysis

Foundations of Modern Probability

Weakest precondition reasoning for expected runtimes of randomized algorithms

Expected runtime analyis by program verification

KoAT: Web interface, binary, Docker image, and examples available at the web site

Semantics of probabilistic programs

Abstraction, Refinement and Proof for Probabilistic Systems

A new proof rule for almostsure termination

Inferring expected runtimes of probabilistic integer programs using expected sizes

Automated termination analysis of polynomial probabilistic programs

Z3: An efficient SMT solver

Bounded expectations: Resource analysis for probabilistic programs

Analyzing innermost runtime complexity of term rewriting by dependency pairs

Reasoning about recursive probabilistic programs

A complete method for the synthesis of linear ranking functions

Markov Decision Processes: Discrete Stochastic Dynamic Programming

Constraint-based linear-relations analysis

Complexity and resource bound analysis of imperative programs using difference constraints

Complexity verification using guided theorem enumeration

TPDB (Termination Problems Data Base

Automatic verification of probabilistic concurrent finite-state programs

Raising expectations: automating expected cost analysis with types

Cost analysis of nondeterministic probabilistic programs

Acknowledgements We thank Carsten Fuhs for discussions on initial ideas.