Inferring Expected Runtimes of Probabilistic Integer Programs Using Expected Sizes

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.


Introduction
There exist several approaches and tools for automatic complexity analysis of nonprobabilistic 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].
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.

Probabilistic Integer Programs
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.
Definition 1 (Bounds). The set of bounds B is the smallest set with PV ∪R ≥0 ⊆ B, and where b 1 , b 2 ∈ B and v ∈ R ≥1 imply b 1 + b 2 , b 1 · b 2 ∈ B and v b1 ∈ B.
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) = v∈Z |v| · pr(v) < ∞. A distribution function d : Let D denote the set of all bounded distribution functions (our implementation supports Bernoulli, uniform, geometric, hypergeometric, and binomial distributions, see [43] for details).
Definition 2 (PIP). (PV, L, GT , 0 ) is a probabilistic integer program with 1. a finite set of program variables PV ⊆ V 2. a finite non-empty set of program locations L 3. a finite non-empty set of general transitions GT . A general transition g is a finite non-empty set of transitions t = ( , p, τ, η, ), consisting of (a) the start and target locations , ∈ L of transition t, (b) the probability p ≥ 0 that transition t is chosen when g is executed, (c) the guard τ ∈ C of t, and (d) the update function η : PV → Z[V] ∪ 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. Then Conf = (L { ⊥ })×(T {t in , t ⊥ })×Σ 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 ⊥ .

Example 3 (PIP). Consider the PIP in
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].

Complexity Bounds
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.

Runtime and Size Bounds
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.
Definition 4 (Runtime and Size Bounds [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 .

Expected Runtime and Size Bounds
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.

Definition 8 (Expected Runtime and Size Bounds).
• 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)).
Example 9 (Expected Runtime and Size Bounds). Our new techniques from 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 g 1 = {t 1 , t 2 } and g 3 = {t 4 }. 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.

Computing Expected Runtime Bounds
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.

Probabilistic Linear Ranking Functions
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.
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.

Inferring Expected Runtime Bounds
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. 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.
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.

Computing Expected Size Bounds
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.

Local Change Bounds and General Result Variable Graph
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 ∈ Σ.
Definition 18 (Expected Local Change Bound). Let SB be a size bound.
The following theorem shows that for any admissible configuration in a state s , CB E (g, , x) is an upper bound on the expected value of |s(x) − s (x)| if s is the next state obtained when applying g in state s to reach location .
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.

Inferring Expected Size Bounds
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.

Inferring Expected Size Bounds for Trivial SCCs 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, , ).
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.2.2 how to compute such expected size bounds for non-trivial SCCs. Hence, by Thm. 23

Inferring Expected Size Bounds for Non-Trivial SCCs
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.
Theorem 25 (Expected Size Bounds for Non-Trivial SCCs). Let (RB E , SB E ) be an expected bound pair, (RB, SB) a (non-probabilistic) bound pair, and let C ⊆ GRV form a non-trivial SCC of the general result variable graph where GT C = {g ∈ GT | (g, , ) ∈ C}. 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.
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].
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. 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).