Automatic Generation of Moment-Based Invariants for Prob-Solvable Loops

TL;DR

This paper introduces a method to automatically generate moment-based invariants for Prob-Solvable loops in probabilistic programs, enabling comprehensive analysis of loop behaviors beyond expected values.

Contribution

It presents a novel combination of symbolic summation and statistical methods to derive higher-order moment invariants for a specific class of probabilistic loops.

Findings

Successfully derived invariants for complex probabilistic loops

Automated computation of higher-order moments demonstrated

Enhanced analysis capabilities for probabilistic program behaviors

Abstract

One of the main challenges in the analysis of probabilistic programs is to compute invariant properties that summarise loop behaviours. Automation of invariant generation is still at its infancy and most of the times targets only expected values of the program variables, which is insufficient to recover the full probabilistic program behaviour. We present a method to automatically generate moment-based invariants of a subclass of probabilistic programs, called Prob-Solvable loops, with polynomial assignments over random variables and parametrised distributions. We combine methods from symbolic summation and statistics to derive invariants as valid properties over higher-order moments, such as expected values or variances, of program variables. We successfully evaluated our work on several examples where full automation for computing higher-order moments and invariants over program variables was not yet possible.

Tables

Table 1. Table 1: Moment-based invariants of Prob-solvable loops, where n 𝑛 n is the loop counter.

Program	Moment	Runtime ($s$)	Computed Moment-Based Invariants
Coupon [18]	1	0.37	$E\left[c\left(n\right)\right]=\left(2^n-1\right)/\left(2^n\right)$
	2	0.40	$E\left[c^2\left(n\right)\right]=\left(2^n-1\right)/\left(2^n\right)$
	3	0.34	$E\left[c^2\left(n\right)\right]=\left(2^n-1\right)/\left(2^n\right)$
Coupon4 [18]	1	0.90	$E\left[c\left(n\right)\right]=\left(4^n-3^3\right)/\left(4^n\right)$
	2	1.1	$E\left[c^2\left(n\right)\right]=\left(4^n-3^3\right)/\left(4^n\right)$
	3	1.3	$E\left[c^3\left(n\right)\right]=\left(4^n-3^3\right)/\left(4^n\right)$
random_walk_1d_cts [18]	1	0.12	$E\left[x\left(n\right)\right]=n/5$
	2	0.45	$E\left[x^2\left(n\right)\right]=n^2/25+22n/75$
	3	1.00	$E\left[x^3\left(n\right)\right]=n^3/125+n^{2}22/125-n21/250$
sum_rnd_series [6]	1	0.31	$E\left[x\left(n\right)\right]=n^2/4+n/4$
	2	2.89	$E\left[x^2\left(n\right)\right]=n^4/16+5n^3/24+3n^2/16+n/24$
	3	17.7	$E\left[x^3\left(n\right)\right]=n^6/64+7n^5/64+13n^4/64+9n^3/64+n^2/32$
product_dep_var [6]	1	0.65	$E\left[p\left(n\right)\right]=n^2/4-n/4$
	2	6.27	$E\left[p^{(}n\right)]=n{}^4/16-n{}^3/8+3n{}^2/16-n/8$
	3	37.5	$E\left[p^3\left(n\right)\right]=n^6/64-3n^5/64+9n^4/64-21n^3/64+15n^2/32-n/4$
random_walk_2d [4, 18]	1	0.07	$E\left[x\left(n\right)\right]=0$
	2	0.26	$E\left[x^2\left(n\right)\right]=n/2$
	3	0.49	$E\left[x^3\left(n\right)\right]=0$
Binomial(”$p$”) [6, 8, 14]	1	0.17	$E\left[x\left(n\right)\right]=np$
	2	0.47	$E\left[x^2\left(n\right)\right]=n^2{p}^2+np\left(1-p\right)$
	3	1.6	$E\left[x^3\left(n\right)\right]=n^3{p}^3-3n^2{p}^3+3n^2{p}^2+2n{p}^3-3n{p}^2+np$
StutteringA – Fig. 1(A)	1	0.44	$E\left[s\left(n\right)\right]=9n/4$
	2	2.2	$E\left[s^2\left(n\right)\right]=81n^2/16+\left(20d^2+27\right)/16n$
	3	8.48	$E\left[s^3\left(n\right)\right]=81d^2{n}^2/16+63d^{2}n/16+729n^3/64+9n^2\left(4d^2-9\right)/32+9n^2\left(4d^2+9\right)/16+567n^2/64+3n\left(-6d^2-21\right)/8+3n\left(6d^2-12\right)/16+243n/32$
StutteringB – Fig. 1(B)	1	0.49	$E\left[s\left(n\right)\right]=9n/4$
	2	2.03	$E\left[s^2\left(n\right)\right]=81n^2/16+347/16n+128/3$
	3	7.43	$E\left[s^3\left(n\right)\right]=729n^3/64+9369n^2/64+1359n/32$=
StutteringC – Fig. 1(C)	1	1.8	$E\left[s\left(n\right)\right]=3n^3/8+3n^2/8-n$
	2	72.5	$E\left[s^2\left(n\right)\right]=9n^6/64+3n^5\left(8d^2+27\right)/160+n^4\left(8d^4+84d^2-90\right)/192+n^3\left(32d^4+216d^2-252\right)/288+n^2\left(8d^4+44d^2+61\right)/64+n\left(80d^4+324d^2-9\right)/1440$
	3	2144	$E\left[s^3\left(n\right)\right]=27n^9/512+27n^8\left(16d^2+39\right)/2560+3n^7\left(824d^4+6444d^2+1242\right)/17920+n^6\left(1900d^4+3996d^2-4365\right)/2560+n^5\left(2004d^4+1704d^2-54\right)/2560+n^4\left(-1900d^4-7056d^2+13446\right)/7680+n^3\left(-6948d^4-12708d^2-6969\right)/7680+n^2\left(-1900d^4-3114d^2-315\right)/3840+n\left(-108d^4-603d^2+288\right)/6720$
StutteringD – Fig. 1(D)	1	1.92	$E\left[s\left(n\right)\right]=3n^3/8+3n^2/8-n$
	2	46.3	$E\left[s^2\left(n\right)\right]=9n^6/64+93n^5/32+1651n^4/96+2849n^3/72+2813n^2/64+5131n/288$
	3	2076	$E\left[s^3\left(n\right)\right]=27n^9/512+1593n^8/512+94587n^7/1792+545971n^6/2560+270117n^5/1280-58585n^4/768-132599n^3/512-536539n^2/3840-771n/140$
StutteringP	1	0.28	$E\left[s\left(n\right)\right]=3np$
	2	1.68	$E\left[s^2\left(n\right)\right]=11n^2{p}^2+3np\left(-2p+1\right)+np\left(-p-1\right)+4np\left(-p+2\right)-1$
	3	6.05	$E\left[s^3\left(n\right)\right]=27n_1^3{p}^3-3n_1^2{p}^3+3n_1^2{p}^2\left(-6p+3\right)+12n_1^2{p}^2\left(-3p+3\right)+12n_1^2{p}^2\left(-2p+3\right)+3n_{1}p\left(4p^2-3p+3\right)+3n_{1}p\left(8p^2-12p+9\right)+n_{1}p\left(p^2-3p\left(-p-1\right)-3p+2\right)/2+2n_{1}p\left(2p^2-6p\left(-p+2\right)-6p+13\right)+6$
Square	1	0.38	$E\left[y\left(n\right)\right]=n^2+n$
	2	2.46	$E\left[y^2\left(n\right)\right]=n^4+6\ast n^3+3\ast n^2-2\ast n$
	3	8.70	$E\left[y^3\left(n\right)\right]=n^6+15\ast n^5+45\ast n^4-15\ast n^3-30\ast n^2+16\ast n$

Equations

$$f(n+r)=R(f(n),f(n+1),\dots ,f(n+r-1),n),\text{with }n\in \mathbb{N},$$

$$f(n+r)=a_{0}f(n)+a_{1}f(n+1)+\dots +a_{r-1}f(n+r-1),\text{with }r,n\in \mathbb{N},$$

$$f(n)=P_1(n){\theta }_1^n+\cdots +P_s(n){\theta }_s^n,$$

$$Mo{m}_k[c,X]=E[(X-c{)}^k]$$

$$M_X(t)=E[e^{tX}],\text{with}t\in \mathbb{R}$$

$$I;\text{while(true)}\{U\},\text{where:}$$

$$x_i:=a_i{x}_i+P_i(x_1,\dots x_{i-1})[p_i]b_i{x}_i+Q_i(x_1,\dots x_{i-1}),$$

$$x_i:=a_i{x}_i+P_i(x_1,\dots x_{i-1}),$$

$$x_i(n+1)=a_i{x}_i(n)+P_i(x_1(n),\dots ,x_{i-1}(n))[p_i]b_i{x}_i(n)+Q_i(x_1(n),\dots ,x_{i-1}(n)).$$

$$E[x(n+r)]=R(E[x(n)],E[x(n+1)],\dots ,E[x(n+r-1)],n)(n\in \mathbb{N}),$$

$$\begin{array}[]{lcl}E[x_{i}(n+1)]&=&E\big{[}p_{i}\cdot f_{1}(x_{i}(n))+(1-p_{i})\cdot f_{2}(x_{i}(n))\big{]}\\ &=&p_{i}\cdot E\big{[}f_{1}(x_{i}(n))\big{]}+(1-p_{i})\cdot E\big{[}f_{2}(x_{i}(n))\big{]}\end{array}$$

$$\begin{array}[]{lcl}E[x_{i}(n+1)]&=&p_{i}\cdot E\big{[}a_{i}x_{i}(n)+P_{i}(x_{1}(n),\dots,x_{i-1}(n))\big{]}\\ &&+(1-p_{i})\cdot E\big{[}b_{i}x_{i}(n)+Q_{i}(x_{1}(n),\dots,x_{i-1}(n))\big{]}.\end{array}$$

$$\begin{array}[]{lcl}E[expr_{1}+expr_{2}]&\rightarrow&E[expr_{1}]+E[expr_{2}]\\ E[expr_{1}\cdot expr_{2}]&\rightarrow&E[expr_{1}]\cdot E[expr_{2}],\ {\small\text{if\ }expr_{1},expr_{2}\text{ are independent}}\\ E[c\cdot expr_{1}]&\rightarrow&c\cdot E[expr_{1}]\\ E[c]&\rightarrow&c\\ E[\mathcal{D}\cdot expr_{1}]&\rightarrow&E[\mathcal{D}]\cdot E[expr_{1}]\end{array}$$

$$\left\{\begin{array}[]{lcl}E[f(n+1)]&=&\frac{3}{4}E[1]+\frac{1}{4}E[0])\\ E[x(n+1)]&=&E\big{[}x(n)+f(n+1)\cdot rand(1-d,1+d)\big{]}\\ E[y(n+1)]&=&E\big{[}y(n)+f(n+1)\cdot rand(2-2d,2+2d)\big{]}\\ E[s(n+1)]&=&E\big{[}x(n+1)+y(n+1)\big{]}\end{array}\right.$$

$$\left\{\begin{array}[]{lcl}E[f(n+1)]&=&\frac{3}{4}\\ E[x(n+1)]&=&E[x(n)]+E[f(n+1)]\cdot E[rand(1-d,1+d)]\\ E[y(n+1)]&=&E[y(n)]+E[f(n+1)]\cdot E[rand(2-2d,2+2d)]\\ E[s(n+1)]&=&E[x(n+1)]+E[y(n+1)]\end{array}\right.$$

$$x_i:=\sum _{k=1}^m{a}_{i,k}{x}_k+c_i[p_i]\sum _{k=1}^m{b}_{i,k}{x}_k+d_i,$$

$$I;\text{while(true)}\{\text{if }t\text{ then }{U}_1\text{ else }{U}_2\},\text{where:}$$

$$\begin{array}[]{l}f:=1[p]0;\\ g:=1[p]0;\\ x:=t(u_{1}f+v_{1}(1-f))+(1-t)(u_{2}g+v_{2}(1-g))\end{array}$$

$$\begin{array}[]{lcl}E[x_{i}(n+1)]&=&p_{i}\cdot E\big{[}a_{i}x_{i}(n)+P_{i}(x_{1}(n),\dots,x_{i-1}(n))\big{]}\\ &&+(1-p_{i})\cdot E\big{[}b_{i}x_{i}(n)+Q_{i}(x_{1}(n),\dots,x_{i-1}(n))\big{]}.\par\end{array}$$

$$Mo{m}_k[0,x_i(n)]=E[x_i(n{)}^k].$$

$$E[x(n+1{)}^2]=E[x(n{)}^2]+2\cdot E[x(n)]\cdot E[f(n+1)]+E[f(n+1{)}^2]\cdot \frac{1}{3}(d^2+3),$$

$$E[y(n+1{)}^2]=E[y(n{)}^2]+4\cdot E[y(n)]\cdot E[f(n+1)]+E[f(n+1{)}^2]\cdot \frac{4}{3}(d^2+3),$$

$$E[s(n+1)^{2}]=E[\big{(}x(n+1)+y(n+1)\big{)}^{2}]=E[x(n+1)^{2}]+2E[(xy)(n+1)]+E[y(n+1)^{2}]$$

$$E[(xy)(n+1)]=E[(xy)(n)]+2\cdot E[x(n)]\cdot E[f(n+1)+E[y(n)]\cdot E[f(n+1)]+2\cdot E[f(n+1{)}^2],$$

$$\begin{array}[]{lclcl}E[f(n+1)]&=&\frac{3}{4}\\ E[x(n+1)]&=&E[x(n)]+E[f(n+1)\cdot{\tt rand(1-d,1+d)}]&=&E[x(n)]+\frac{3}{4}\\ E[y(n+1)]&=&E[y(n)]+E[f(n+1)\cdot{\tt rand(2-2d,2+2d)}]&=&E[x(n)]+2\cdot\frac{3}{4}\\ E[s(n+1)]&=&E[x(n+1)]+E[y(n+1)]\end{array}$$

$$\begin{array}[]{lclclcl}E[f(n)]&=&\frac{3}{4}&&E[s(n)]&=&\frac{9}{4}n\\ E[x(n)]&=&\frac{3}{4}n-1&&E[y(n)]&=&\frac{3}{2}n+1\end{array}$$

$$\begin{array}[]{lclclcl}E[f(n)^{2}]&=&\frac{3}{4}&&E[s(n)^{2}]&=&\frac{81}{16}n^{2}+\frac{20d^{2}+27}{16}n\\ E[x(n)^{2}]&=&\frac{9}{16}n^{2}+\frac{4d^{2}-21}{16}n+1&&E[y(n)^{2}]&=&\frac{9}{4}n^{2}+\frac{4d^{2}+15}{4}n+1\\ E[(xy)(n)]&=&\frac{9}{8}n^{2}-\frac{3}{8}n-1\end{array}$$

$$x_k^{{\alpha }_k}\cdot x_{k-1}^{{\alpha }_{k-1}}\dots x_1^{{\alpha }_1}\stackrel{\sigma }{\to }{\omega }^k\cdot {\alpha }_k+{\omega }^{k-1}\cdot {\alpha }_{k-1}\cdots +{\alpha }_1,$$

$$x_{i}^{\alpha_{i}}:=\big{(}c_{i}x_{i}+P_{i}(x_{1},\cdots x_{i-1})\big{)}^{\alpha_{i}}$$

$$E[M(n+1)]$$

Full text

1 Introduction

Probabilistic programs (PPs), originally employed in cryptographic/privacy protocols and randomised algorithms, are now gaining momentum due to the several emerging applications in the areas of machine learning and AI [10]. One of the main problems that arise from introducing randomness into the program is that we can no longer view variables as single values; we must think about them as distributions. Existing approaches, see e.g. [21, 14, 4] usually take into consideration only expected values, or upper and lower bounds over program variables. As argued by [23], such information is however insufficient to characterize the full value distributions of variables; (co-)variances and other higher-order moments of variables are also needed. Consider for example the PPs of Fig. 1(A) and Fig. 1(B): the expected value of variable $s$ at each loop iteration is the same in both PPs, while the variance of the value distribution of $s$ differs in general (a similar behaviour is also exploited by Fig. 1(C)-(D)). Thus, Fig. 1(A) and Fig. 1(B) do not have the same invariants over higher-order moments; yet, current approaches would fail identifying such differences and only compute expected values of variables.

One of the main challenges in analysing PPs and computing their higher-order moments comes with the presence of loops and the burden of computing so-called quantitative invariants [14]. Quantitative invariants are properties that are true before and after each loop iteration. Weakest pre-expectations [21, 14] can be used to compute quantitative invariants. This approach, supported for example in  Prinsys [11], consists in annotating a loop with a template invariant and then solve constraints over the unknown coefficients of the template. Other methods [2, 18] use martingales that are expressions over program variables whose expectations remain invariant. The aforementioned approaches are however not fully automatic since they require user guidance for providing templates and hints. In addition, they are limited to invariants over only expected values: with the exception of [18], they do not compute higher-order moments describing the distribution generated by the PP (see Section 6 for more details).

In this paper we introduce a fully automated approach to compute invariant properties over higher-order moments of so-called Prob-solvable loops, to stand for probabilistic P-solvable loops. Prob-solvable loops are PPs that extend the imperative P-solvable loops described in [17] with probabilistic assignments over random variables and parametrised distributions. As such, variable updates are expressed by random polynomial, and not only affine, updates (see Section 3). Each program in Fig. 1 is Prob-solvable; moreover, Fig. 1(C)-(D) involve nonlinear updates over $s$.

Our work uses statistical properties to eliminate probabilistic choices and turn random updates into recurrence relations over higher-order moments of program variables. We show that higher-order moments of Prob-solvable loops can be described by C-finite recurrences (Theorem 4.1). We further solve such recurrences to derive moment-based invariants of Prob-solvable loops (Section 4). A moment-based invariant is a property that holds at arbitrary loop iterations (hence, invariants), expressing closed form solutions of higher-order moments of program variables. To the best of our knowledge, no other method is able to derive higher-order moments of PPs in a fully automated approach. Our work hence allows to replace, for example, the required human guidance of [11, 19] for Prob-solvable loops. Unlike existing works, we support PPs with parametrised distributions (e.g., in Fig. 1(A)): instead of taking concrete instances of a given parametrised distribution, we automatically infer invariants of the entire class of PPs characterised by the considered parametrised distribution.

Our approach is both sound and terminating: given a Prob-solvable loops and an integer $k\geq 1$, we automatically infer the moment-based invariants over the $k$th moments of our input loop (see Section 4). Unlike the approach of [17] for deriving polynomial invariants of non-probabilistic (P-solvable) loops, our work only computes closed form expressions over higher-order moments and does not employ Gröbner basis computation to eliminate loop counters from the derived closed forms. As such, our moment-based invariants are not restrictive to polynomial properties but are linear combinations of polynomial expressions and exponential sequences over the loop counter. Moreover, Prob-solvable are more expressive than P-solvable loops as they are not restricted to deterministic updates but allow random assignments over variables.

Contributions.

Our main contributions are: (1) we introduce the class of Prob-solvable loops with probabilistic assignments over random variables and distributions (Section 3); (2) we show that Prob-solvable loops can be modelled as C-finite recurrences over higher-order moments of variables (Theorem 4.1); (3) we provide a fully automated approach that derives moment-based invariants over arbitrary higher-order moments of Prob-solvable loops (Algorithm 1); (4) we implemented our work as an extension of the Aligator package [12] and evaluated over several challenging PPs (Section 5).

2 Preliminaries

We recall basic mathematical properties about recurrences and higher-order moments of variable values – for more details see [16, 20]. Throughout this paper, let $\mathbb{N},\mathbb{Z},\mathbb{R}$ denote the set of natural, integer and real numbers. We reserve capital letters to denote abstract random variables, e.g. $X,Y,\ldots$, and use small letters to denote program variables, e.g. $x,y,\ldots$, all possibly with indices.

2.1 C-Finite Recurrences

While sequences and recurrences are defined over arbitrary fields of characteristic zero, in our work we only focus over sequences/recurrences over $\mathbb{R}$.

Definition 1 (Sequence)

A univariate sequence in $\mathbb{R}$ is a function $f:\mathbb{Z}\rightarrow\mathbb{R}$. A recurrence for a sequence $f(n)$ is

[TABLE]

for some function $R:\mathbb{R}^{r+1}\rightarrow\mathbb{R}$, where $r\in\mathbb{N}$ is called the order of the recurrence.

For simplicity, we denote by $f(n)$ both the recurrence of $f(n)$ as well as the recurrence equation $f(n)=0$. When solving the recurrence $f(n)$, one is interested in computing a closed form solution of $f(n)$, expressing the value of $f(n)$ as a function of $n$ for any $n\in\mathbb{N}$. In our work we only consider the class of linear recurrences with constant coefficients, also called C-finite recurrences.

Definition 2 (C-finite recurrences)

A C-finite recurrence $f(n)$ satisfies the linear homogeneous recurrence with constant coefficients:

[TABLE]

where $r$ is the order of the recurrence, and $a_{0},\dots,a_{r-1}\in\mathbb{R}$ are constants with $a_{0}\neq 0$.

An example of a C-finite recurrence is the recurrence of Fibonacci numbers satisfying the recurrence $f(n+2)=f(n+1)+f(n)$, with initial values $f(0)=0$ and $f(1)=1$. Unlike arbitrary recurrences, closed forms of C-finite recurrences $f(n)$ always exist [16] and are defined as:

[TABLE]

where $\theta_{1},\ldots,\theta_{s}\in{\mathbb{R}}$ are the distinct roots of the characteristic polynomial of $f(n)$ and $P_{i}(n)$ are polynomials in $n$. Closed forms of C-finite recurrences are called C-finite expressions. We note that, while the C-finite recurrence (1) is homogeneous, inhomogeneous C-finite recurrences can always be translated into homogeneous ones, as the inhomogeneous part of a C-finite recurrence is a C-finite expression.

In our work, we focus on the analysis of Prob-solvable loops and consider loop variables $x$ as sequences $x(n)$, where $n\in\mathbb{N}$ denotes the loop iteration counter. Thus, $x(n)$ gives the value of the program variable $x$ at iteration $n$.

2.2 Expected Values and Moments of Random Variables

Here we introduce the relevant notions from statistics that our work relies upon.

Definition 3 (Probability space)

A probability space is a triple $(\Omega,F,P)$ consisting of a sample space $\Omega$ denoting the set of outcomes, where $\Omega\not=\emptyset$, a $\sigma$-algebra $F$ with $F\subset 2^{\Omega}$, denoting a set of events, a probability measure $P:F\rightarrow[0,1]$ s.t. $P(\Omega)=1$.

We now define random variables and their higher-order moments.

Definition 4 (Random variable)

A random variable $X:\Omega\rightarrow\mathbb{R}$ is a measurable function from a set $\Omega$ of possible outcomes to $\mathbb{R}$.

In the context of our Prob-solvable loops, for each loop variable $x$, we consider elements of its corresponding sequence $x(n)$ to be random variables. When working with a random variable $X$, one is in general interested in expected values and other moments of $X$.

Definition 5 (Expected value)

An expected value of a random variable $X$ defined on a probability space $(\Omega,F,P)$ is the Lebesgue integral: $E[X]=\int_{\Omega}X\cdot dP.$ In the special case when $\Omega$ is discrete, that is the outcomes are ${X_{1},\dots X_{N}}$ with corresponding probabilities ${p_{1},\dots p_{N}}$, we have $E[X]=\sum_{i=1}^{N}X_{i}\cdot p_{i}.$ The expected value of $X$ is often also referred to as the mean or $\mu$ of $X$.

For program variables $x$ of Prob-solvable loops, our work computes the expected values of the corresponding sequences $x(n)$ but also higher-order and mixed moments.

Definition 6 (Higher-Order Moments)

Let $X$ be a random variable, $c\in\mathbb{R}$ and $k\in\mathbb{N}$. We write $Mom_{k}[c,X]$ to denote the $k$th moment about $c$ of $X$, which is defined as:

[TABLE]

In this paper we will be almost solely interested in moments about [math] (called raw moments) and about the mean $E[X]$ (called central moments). We note though that we can move to moments of $X$ with different centers using Proposition 1.

Proposition 1 (Transformation of center)

Let $X$ be a random variable, $c,d\in\mathbb{R}$ and $k\in\mathbb{N}$. The $k$th moment about $d$ of $X$, can be calculated from moments about $c$ of $X$ by: $\displaystyle E\left[(X-d)^{k}\right]=\sum_{i=0}^{k}{k\choose i}E\left[(X-c)^{i}\right](c-d)^{k-i}.$

Similarly to higher-order moments, we also consider mixed moments, that is $E[X\cdot Y]$, where $X$ and $Y$ are random variables. For arbitrary random variables $X$ and $Y$, we have the following basic properties about their expected values and other moments:

• •

$E[c]=c$ for a constant $c\in\mathbb{R}$,

• •

expected value is linear, $E[X+Y]=E[X]+E[Y]$ and $E[c\cdot X]=c\cdot E[X]$,

• •

expected value is not multiplicative, in general $E[X\cdot Y]\not=E[X]\cdot E[Y]$

• •

expected value is multiplicative for independent random variables.

As a consequence of the above, expected values of monomials over arbitrary random variables, e.g. $E[X\cdot Y^{2}]$, cannot be in general further simplified.

The moments of a random variable $X$ with bounded support fully characterise its value distribution. While computing all moments of $X$ is generally very hard, knowing only a few moments of $X$ gives useful information about its value distributions. The most common moments are variance, covariance, skewness, as defined below.

Definition 7 (Common moments)

Variance measures how spread the distribution is and is defined as the second central moment: $Var[X]=Mom_{2}[E[X],X]$.

Covariance is a mixed moment measuring variability of two distributions and is defined as: $Cov[X,Y]=E[(X-E[X])\cdot(Y-E[Y])]$.

Skewness measures asymmetry of the distribution and is defined as the normalised third central moment: $Skew[X]=\frac{Mom_{3}[E[X],X]}{(Var[X])^{3/2}}$.

Basic results about variance and covariance state: $Cov[X,X]=Var[X]$, $Var[X]=E[X^{2}]-(E[X]^{2})$, and $Cov[X,Y]=E[X\cdot Y]-E[X]\cdot E[Y]$.

Definition 8 (Moment-Generating Function (MGF))

A moment generating function of a random variable $X$ is given by:

[TABLE]

whenever this expectation exists.

Moment-generating functions, as the name suggests, can be used to compute higher-order moments of a random variable $X$. If we take the $k$th derivative of the moment-generating function of $X$, evaluated at [math], we get the $k$th moment about [math] of $X$, that is $Mom_{k}[0,X]$111due to the series expansion $e^{tX}=1+tE[X]+\frac{t^{2}E[X^{2}]}{2!}+\frac{t^{3}E[X^{3}]}{3!}+\dots$ and derivative w.r.t. $t$. For many standard distributions, including Bernoulli, uniform and normal distributions, the moment-generating function exists and gives us a way to compute the moments for random variables drawing from these distributions. Thanks to these properties, we can use common distributions in our Prob-solvable programs.

3 Programming Model: Prob-solvable Programs

We now introduce our programming model of Prob-solvable programs, to stand for probabilistic P-solvable programs. P-solvable programs [17] are non-deterministic loops whose behaviour can be expressed by a system of C-finite recurrences over program variables. Prob-solvable programs extend P-solvable programs by allowing probabilistic assignments over random variables and distributions.

Prob-solvable Loops.

Let $m\in\mathbb{N}$ and $x_{1},\ldots x_{m}$ denote real-valued program variables. We define Prob-solvable loops with $x_{1},\ldots x_{m}$ variables as programs of the form:

[TABLE]

• •

$I$ is a sequence of initial assignments over $x_{1},\ldots,x_{m}$. That is, $I$ is an assignments sequence $x_{1}:=c_{1};x_{2}:=c_{2};\dots x_{m}:=c_{m}$, with $c_{i}\in\mathbb{R}$ representing a number drawn from a known distribution 222a known distribution is a distribution with known and computable moments - in particular, $c_{i}$ can be a real constant.

• •

$U$ is the loop body and is a sequence of $m$ random updates, each of the form:

[TABLE]

or, in case of a deterministic assignment,

[TABLE]

where $a_{i},b_{i}\in\mathbb{R}$ are constants and $P_{i},Q_{i}\in\mathbb{R}[x_{1},\ldots,x_{i-1}]$ are polynomials over program variables $x_{1},\ldots,x_{i-1}$. Further, $p_{i}\in[0,1]$ in (6) is the probability of updating $x_{i}$ to $a_{i}x_{i}+P_{i}(x_{1},\dots x_{i-1})$, whereas the probability to update $x_{i}$ to $b_{i}x_{i}+Q_{i}(x_{1},\dots x_{i-1})$ in (6) is $1-p_{i}$.

The coefficients $a_{i}$, $b_{i}$ and the coefficients of $P_{i}$ and $Q_{i}$ in the variable assignments (6)-(7) of Prob-solvable loops can be drawn from a random distribution as long as the moments of this distribution are known and are independent from program variables $x_{1},\ldots,x_{m}$. Hence, the variable updates of Prob-solvable loop can involve coefficients drawn from Bernoulli, uniform, normal, and other distributions. Moreover, Prob-solvable support parametrised distributions, for example one may have the random distribution ${\tt rand(d_{1},d_{2})}$ with arbitrary $d_{1},d_{2}\in\mathbb{R}$ symbolic constants. Similarly, rather than only considering concrete numeric values of $p_{i}$, the probabilities $p_{i}$ in the probabilistic updates (6) of Prob-solvable loops can also be symbolic constants.

Example 1

The programs in Fig. 1 are Prob-solvable, using uniform distributions given by ${\tt rand()}$. Fig. 1(D) also uses normal distribution given by ${\tt gauss()}$. Note that while the random distributions of Fig. 1(B,D) are defined in terms of concrete constants, Fig. 1(A,C) have a parametrised random distribution, defined in terms of $d\in\mathbb{R}$.

Prob-solvable Loops and Moment-Based

Recurrences.

Let us now consider a Prob-solvable program with $n\in\mathbb{N}$ denoting the loop iteration counter. We show that variable updates of Prob-solvable programs yield special recurrences in $n$, called moment-based recurrences. For this, we consider program variables $x_{1},\ldots,x_{m}$ as sequences $x_{1}(n),\ldots,x_{m}(n)$ allowing us to precisely describe relations between values of $x_{i}$ at different loop iterations. Using this notation, probabilistic updates (6) over $x_{i}$ turn $x_{i}(n)$ into a random variable, yielding the relation (similarly, for deterministic updates (7)):

[TABLE]

The relation above could be treated as a recurrence equation over random variables $x_{i}(n)$ provided the probabilistic behaviour depending on $p$ is encoded (as an extension) into a recurrence equation. To analyse such probabilistic updates of Prob-solvable loops, for each random variable $x_{i}(n)$ we consider their expected values $E[x_{i}(n)]$ and create new recurrence variables from expected values of monomials over original program variables (e.g. a new variable $E[x_{i}\cdot x_{j}]$). We refer to these new recurrence variables as E-variables. We note that any program variable yields an E-variable, but not every E-variable corresponds to one single program variable as E-variables are expected values of monomials over program variables. We now formulate recurrence equations over E-variables rather than over program variables, yielding moment-based recurrences.

Definition 9 (Moment-Based Recurrences)

Let $x(n)$ be a sequence of random variables. A moment-based recurrence for $x$ is a recurrence over E-variable $E[x]$:

[TABLE]

for some function $R:\mathbb{R}^{r+1}\rightarrow\mathbb{R}$, where $r\in\mathbb{N}$ is the order of the moment-based recurrence.

Similarly to [21], note that variable updates $x_{i}:=f_{1}(x_{i})\>[p_{i}]\>f_{2}(x_{i})$ yield the relation:

[TABLE]

Thanks to this relation, probabilistic updates (6) are rewritten into the moment-based recurrence equation:

[TABLE]

In particular, we have $E[x_{i}(n+1)]=p_{i}\cdot E[a_{i}x_{i}(n)+P_{i}(x_{1}(n),\dots,x_{i-1}(n))]$ for the deterministic assignments of (7) (that is, $p_{i}=1$ in (7)).

By using properties of expected values of expressions $expr_{1},expr_{2}$ over random variables, we obtain the following simplification rules:

[TABLE]

where $c\in\mathbb{R}$ is a constant and $\mathcal{D}$ is a known independent distribution.

Example 2

The moment-based recurrences of the Prob-solvable loop of Fig. 1(A) are:

[TABLE]

By using the simplification rules (10) on the above recurrences, we obtain the following simplified moment-based recurrences of Fig. 1(A):

[TABLE]

In Section 4 we show that Prob-solvable loops can further be rewritten into a system of C-finite recurrences over E-variables.

Prob-solvable Loops and Mutually Dependent

Updates.

Consider PP loops with mutually dependent affine updates:

[TABLE]

where $a_{i,k},b_{i,k},c_{i},d_{i}\in\mathbb{R}$ are constants. While such assignments are not directly captured by updates (6) of Prob-solvable loops, this is not a restriction of our work. Variable updates given by (12) yield mutually dependent C-finite recurrences over E-variables. Using methods from [16], this coupled system of C-finite recurrences can be rewritten into an equivalent system of independent C-finite recurrences over E-variables, yielding an independent system of moment-based recurrences over which our invariant generation algorithm from Section 4 can be applied. Hence probabilistic loops with affine updates are special cases of Prob-solvable loops.

Multi-Path Prob-solvable Loops.

While (5) defines Prob-solvable programs as single-path loops, the following class of multi-path loops can naturally be modeled by Prob-solvable programs:

[TABLE]

$I$ is as in (5), $t$ is a boolean-valued random variable, and $U_{1}$ and $U_{2}$ are respectively sequences of deterministic updates $x_{i}:=a_{i}x_{i}+P_{i}(x_{1},\dots x_{i-1})$ and $x_{i}:=b_{i}x_{i}+Q_{i}(x_{1},\dots x_{i-1})$ as in (7). PPs (13) can be rewritten to equivalent Prob-solvable loops, as follows. A pair of updates $x:=u_{1}[p]v_{1}$ from $U_{1}$ and $x:=u_{2}[p]v_{2}$ from $U_{2}$ is rewritten by the following sequence of updates:

[TABLE]

with $f,g$ fresh program variables. The resulting program is Prob-solvable and we can thus compute moment-based invariants of multi-path loops as in (13). The programs Coupon, Random_Walk_2D of Table 1 are Prob-solvable loops corresponding to such multi-path loops.

4 Moment-Based Invariants of Prob-solvable Loops

Thanks to probabilistic updates, the values of program variables of Prob-solvable loops after specific number of loop iterations are not a priory determined. The value distributions $x_{i}(n)$ of program variables $x_{i}$ are therefore random variables. When analysing Prob-solvable loops, and in general probabilistic programs, one is therefore required to capture relevant properties over expected values and higher moments of the variables in order to precisely summarise the value distribution of program variables.

Moment-Based Invariants.

We are interested in automatically generating so-called moment-based invariants of Prob-solvable loops. Moment-based invariants are properties over expected values and higher moments of program variables such that these properties hold at arbitrary loop iterations (and hence are invariants).

Automated Generation of Moment-Based

Invariants of Prob-solvable Loops.

Our method for generating moment-based invariants of Prob-solvable loops is summarized in Algorithm 1. Algorithm 1 takes as input a Prob-solvable loop $\mathcal{P}$ and a natural number $k\geq 1$ and returns moment-based invariants over the $k$th moments of the program variables $\{x_{1},\ldots,x_{m}\}$. We denote by $n$ the loop counter of $\mathcal{P}$.

Theorem 4.1

Higher-order moments of variables in Prob-solvable loops can be modeled by C-finite recurrences over E-variables.

Proof (Sketch)

We want to show that $E[x_{i}^{\alpha_{i}}]$ can be expressed using recurrence equation. The idea is to express $x_{i}^{\alpha_{i}}(n+1)$ in terms of value of $x_{i}$ at $n$-th iteration. Value of $x_{i}(n+1)$ is $a_{i}x_{i}(n)+P_{i}(x_{1}(n+1),\dots x_{i-1}(n+1)$ with probability $p$ and $b_{i}x_{i}(n)+Q_{i}(x_{1}(n+1),\dots x_{i-1}(n+1)$ with probability $(1-p)$. From here we can derive that E[x_{i}^{\alpha_{i}}(n+1)]=E[p_{i}\cdot\big{(}a_{i}x_{i}+P_{i}(x_{1},\dots x_{i-1})\big{)}^{\alpha_{i}}+(1-p_{i})\cdot\big{(}b_{i}x_{i}+Q_{i}(x_{1},\dots x_{i-1})\big{)}^{\alpha_{i}}]. For arbitrary monomial $M=\prod x_{i}^{\alpha_{i}}(n+1)$ we can express $E[M]$ by substituting each $x_{i}^{\alpha_{i}}(n+1)$ as above. This process is captured by line 8 of Alg. 1. The new equations can be further simplified using properties of expected values and the simplification rules (10) to give recurrence equations over E-variables.

We now describe Algorithm 1. Our algorithm first rewrites $\mathcal{P}$ into a set $MBRecs$ of moment-based recurrences, as described in Section 3. That is, program variables $x_{i}$ are turned into random variables $x_{i}(n)$ and variable updates over $x_{i}$ become moment-based recurrences over E-variables by using the relation of (8) (lines 1-2) of Alg. 1).

The algorithm next proceeds with computing the moment-based recurrences of the $k$th moments of $x_{1},\ldots,x_{m}$. Recall that the $k$th moment of $x_{i}$ is given by:

[TABLE]

Hence, the set $S$ of monomials yielding E-variables for which moment-based recurrences need to be solved is initialized to $\{x_{1}^{k},\ldots,x_{m}^{k}\}$ (line 3 of Alg. 1). Note that by considering the resulting E-variables $E[x_{i}^{k}]$ and solving the moment-based recurrences of $E[x_{i}^{k}]$, we derive closed forms of the $k$th moments of $\{x_{1},\ldots,x_{m}\}$ (line 23 of Alg. 1). To this end, Algorithm 1 recursively computes the moment-based recurrences of every E-variable arising from the moment-based recurrences of $E[x_{i}^{k}]$ (lines 5-21 of Alg. 1), thus ultimately computing closed forms for $E[x_{i}^{k}]$. One can then use transformations described in Proposition 1 to compute closed forms for other moments, such as variance and covariance. In more detail,

• •

for each monomial $M=\prod x_{j}^{\alpha_{j}}$ from $S$, we substitute $x_{i}^{\alpha_{i}}$ in $M$ by its probabilistic behaviour. That is, the update of $x_{i}$ in the Prob-solvable loop $\mathcal{P}$ is rewritten, according to (8), into the sum of its two probabilistic updates, weighted by their respective probabilities (lines 6-8 of Alg. 1). Rewriting in line 8 of Alg. 1 represents the most non-trivial step in our algorithm, combining non-deterministic nature of our program with polynomial properties. The resulting polynomial $M^{\prime}$ from $M$ is then reordered to be expressed as a sum of new monomials $N_{j}$ (line 13 of Alg. 1); such a sum always exists as $M^{\prime}$ involves only addition and multiplication over $x_{1},\ldots,x_{m}$ (recall that $P_{i}$ and $Q_{i}$ are polynomials over $x_{1},\ldots,x_{m}$).

• •

By applying the simplification rules(10) of E-variables over the moment-based recurrence of $E[\sum N_{j}]$, the recurrence of $E[M(n+1)]$ is obtained and added to the set $MBRecs$. Here, $M(n+1)$ denotes $\prod_{i=1}^{m}x_{i}(n+1)^{\alpha_{i}}$. As the recurrence of $E[M(n+1)]$ depends on $E[N_{j}]$, moment-based recurrences of $E[N_{j}]$ need also be computed and hence $S$ is enlarged by $N_{j}$ (lines 14-20 of Alg. 1).

As a result, the set $MBRecs$ of moment-based recurrences of E-variables corresponding to $S$ are obtained. These recurrences are C-finite expressions over E-variables (see correctness argument of Theorem 4.3) and hence their closed form solutions exist. In particular, the closed forms $CF_{i}(k,n)$ of $E[x_{i}(n)^{k}]$ is derived, turning $E[x_{i}(n)^{k}]-CF_{i}(k,n)=0$ into a inductive property that holds at arbitrary loop iterations and is hence a moment-based invariant of $\mathcal{P}$ over the $k$th moment of $x_{i}$ (line 23 of Alg. 1).

Theorem 4.2 (Soundness)

Consider a Prob-solvable loop $\mathcal{P}$ with program variables $x_{1},\ldots,x_{m}$ and let $k$ be a non-negative integer with $k\geq 1$. Algorithm 1 generates moment-based invariants of $\mathcal{P}$ over the $k$th moments of $x_{1},\ldots,x_{m}$.

Note when $k=1$, Algorithm 1 computes the moment-based invariants as invariant relations over the closed form solutions of expected values of $x_{1},\ldots,x_{m}$. In this case, our moment-based invariants are quantitative invariants as in [14].

Example 3

We illustrate Algorithm 1 for computing the second moments (i.e. $k=2$) of the Prob-solvable loop of Fig. 1(A). Our algorithm initializes $MBRecs=\{E[f(n+1)],E[x(n+1)],E[y(n+1)],E[s(n+1)]\}$ and $S=\{f^{2},x^{2},y^{2},s^{2}\}$.

We next (arbitrarily) choose $M$ to be the monomial $f^{2}$ from $S$. Thus, $S=\{x^{2},y^{2},s^{2}\}$. Using the probabilistic update of $f$, we replace $f^{2}$ by $\frac{3}{4}\cdot 1^{2}+(1-\frac{3}{4})\cdot^{2}$, that is by $\frac{3}{4}$. As a result, $MBRecs=MBRecs\cup\{E[f(n+1)^{2}]=\frac{3}{4}\}$ and $S$ remains unchanged.

We next choose $M$ to be $x^{2}$ and set $S=\{y^{2},s^{2}\}$. We replace $x^{2}$ by its randomised behaviour, yielding E[M(n+1)]=E[x(n+1)^{2}]=E[\big{(}x(n)+f(n+1)\cdot{\tt rand(1-d,1+d)}\big{)}^{2}]. By the simplification rules (10) over E-variables, we obtain:

[TABLE]

as $f(n+1)$ is independent from $x(n)$ and $E[{\tt rand(1-d,1+d)}^{2}]=\frac{1}{3}(d^{2}+3).$ We add the recurrence (15) to $MBRecs$ and keep $S$ unchanged as the E-variables $E[x(n)],E[f(n+1)],E[f(n+1)^{2}]$ have their recurrences already in $MBRecs$.

We next set $M$ to $y^{2}$ and change $S=\{s^{2}\}$. Similarly to $E[x(n+1)^{2}]$, we get:

[TABLE]

by using that $f(n+1)$ is independent from $y(n)$ and $E[{\tt rand(2-2d,2+2d)}^{2}]=\frac{4}{3}(d^{2}+3)$. We add the recurrence (16) to $MBRecs$ and keep $S$ unchanged.

We set $M$ to $s^{2}$, yielding $S=\emptyset$. We extend $MBRecs$ with the recurrence:

[TABLE]

and add $xy$ to $S$. We therefore consider $M$ to be $xy$ and set $S=\emptyset$. We obtain:

[TABLE]

by using that $E[{\tt rand(1-d,1+d)}]=1$ and $E[{\tt rand(2-2d,2+2d)}]=2$. We add the recurrence of $E[(xy)(n+1)]$ to $MBRecs$ and keep $S=\emptyset$.

As a result, we proceed to solve the moment-based recurrences of $MBRecs$. We focus first on the recurrences over expected values:

[TABLE]

Note that the above recurrences are C-finite recurrences over E-variables. For computing closed forms, we respectively substitute $E[f(n+1)$ by its closed form in $E[y(n+1)]$ and $E[x(n+1)]$, yielding closed forms for $E[y(n+1)]$ and $E[x(n+1)]$, and hence for $E[s(n+1)]$. By also using the the initial values of Fig. 1, we derive the closed forms:

[TABLE]

We next similarly derive the closed forms for higher-order and mixed moments:

[TABLE]

yielding hence the moment-based invariants over the second moments of variables of Fig. 1. Using Proposition 1 and Definition 7, we derive the variance of $s(n)$ as $Var(s(n))=\frac{20d^{2}+27}{16}n.$ ∎

Let us finally note that the termination of Algorithm 1 depends whether for every monomial $M$ (from the set $S$, line 5 of Alg. 1) the moment-based recurrence equation over the corresponding E-variable $E[M(n+1)]$ can be computed as C-finite expression over E-variables. To prove this, one can use transfinite induction over monomials . That is, we order monomials and show that the recurrence of $E[M(n+1)]$ depends only on smaller monomials, for which we can compute C-finite closed form expressions. Thus we have an inhomogeneous C-finite recurrence relation for $E[M(n+1)]$, yielding a C-finite closed form expression.

Theorem 4.3 (Termination)

*For any non-negative integer $k$ with $k\geq 1$ and any Prob-solvable loop $\mathcal{P}$ with program variables $x_{1},\ldots,x_{m}$, Algorithm 1 terminates.

Moreover, Algorithm 1 terminates in at most $\mathcal{O}(k^{m}\cdot d_{m}^{m-1}\cdot d_{m-1}^{m-2}\ldots\cdot d_{2}^{1})$ steps, where $d_{i}=max\{deg(P_{i}),deg(Q_{i}),1\}$ with $deg(P_{i}),deg(Q_{i})$ denoting the degree of polynomials $P_{i}$ and $Q_{i}$ of the variable updates (6).*

Proof

We associate every monomial with an ordinal number as follows:

[TABLE]

and order monomials $M,N$ such that $M>N$ if $\sigma(M)>\sigma(N)$. Algorithm 1 terminates if for every monomial $M$ (from the set $S$, line 5 of Alg. 1) the moment-based recurrence equation over the corresponding E-variable $E[M(n+1)]$ can be computed as C-finite expression over E-variables. We will show that this is indeed the case by transfinite induction over monomials.

Let $M=\prod_{k=1}^{K}x_{k}^{\alpha_{k}}$ be a monomial and assume that every smaller monomial has a closed form solution in form of a C-finite expression.

Let

[TABLE]

be the updates of our variables after removing the probabilistic choice as in line 6 of the algorithm. Then recurrence for $M$ is

[TABLE]

for some $J$, constants $b_{i}$ and monomials $N_{1},\dots,N_{J}$ all different than $M$. By Lemma 1, we have an inhomogeneous C-finite recurrence relation $E[M(n+1)]=E[M(n)]+\gamma$, for some C-finite expression $\gamma$. Hence, the closed form of $E[M(n+1)]$ exists and is a C-finite expression. ∎

We finally prove our auxiliary lemma.

Lemma 1

$M>N_{j}$* for all $j\leq J$ in (4).*

Proof

Let $M=\prod_{k=1}^{K}x_{k}^{\alpha_{k}}$ and have $N_{j}=\prod_{k=1}^{K}x_{k}^{\beta_{k}}$ coming from

[TABLE]

Assume $M\leq N_{j}$, i.e. $\omega^{K}\cdot\alpha_{K}+\dots+\alpha_{1}\leq\omega^{K}\cdot\beta_{K}+\dots+\beta_{1}$, so we have $\alpha_{K}\leq\beta_{K}$. Note that in (19) $x_{K}$ only appears in factor $c_{K}x_{K}+P_{K}(x_{1},\dots x_{K-1})$. Considering the multiplicity, we get at most $\alpha_{K}$th power of $x_{K}$, hence $\alpha_{K}\geq\beta_{K}$. Thus $\alpha_{K}=\beta_{K}$

So for $M\leq N_{j}$ we need $N_{j}$ from (c_{K}x_{K})^{\alpha_{K}}\cdot\prod_{i=1}^{K-1}\big{(}c_{i}x_{i}+P_{i}(x_{1},\cdots x_{i-1})\big{)}^{\alpha_{i}}.

Proceeding similarly for $x_{K-1},x_{K-2},\cdots$ we get that for each $k\leq K$ we have $\alpha_{k}=\beta_{k}$, which contradicts the assumption, thus $M>N_{j}$ as needed.

Regarding the termination, let’s look at what monomials can possibly be added to $S$. Let $M=\prod x_{i}^{\alpha_{i}}\in S$. Based on the algorithm and above it is clear that in case $i=m$ we have $\alpha_{m}\leq k$. For any $i<m$ the maximum value of $\alpha_{i}$ is $\alpha_{i+1}\cdot d_{i+1}$. Hence we have $\alpha_{i}\leq k\cdot\prod_{j=i+1}^{m}d_{j}$. thus we can count all possible monomials, hence the upper bound on the algorithm time complexity, as product of theses upper bounds. This yields $k^{m}\cdot d_{m}^{m-1}\cdot d_{m-1}^{m-2}\ldots\cdot d_{2}^{1}$ as claimed. ∎

5 Implementation and Experiments

We implemented our work in the Julia language, using Aligator[12] for handling and solving recurrences. We evaluated our work on several challenging probabilistic programs with parametrised distributions, symbolic probabilities and/or both discrete and continuous random variables. All our experiments were run on MacBook Pro 2017 with 2.3 GHz Intel Core i5 and 8GB RAM. Our implementation and benchmarks are available at: github.com/miroslav21/aligator.

Benchmarks.

We evaluated our work on 13 probabilistic programs, as follows. We used 7 programs from works [6, 14, 4, 8, 18] on invariant generation. These examples are given in lines 1-7 of Table 1; we note though that Binomial(”$p$”) represents our generalisation of a binomial distribution example taken from [6, 8, 14] to a probabilistic program with parametrised probability $p$. We further crafted 6 examples of our own, illustrating the distinctive features of our work. These examples are listed in lines 8-13 of Table 1: lines 8-11 correspond to the examples of Fig. 1; line 12 of Table 1 shows a variation of Fig. 1, with a parametrized distribution $p$; line 13 corresponds to a non-linear Prob-solvable loop computing squares. All our benchmarks are also available at the aforementioned url.

5.1 Coupon

Probabilistic model of Coupon Collector’s program for two coupons, taken from [18].

⬇

f := 0

c := 0

d := 0

while true:

f := 1 [1/2] 0

c := 1 - f + c*f

d := d + f - d*f

5.2 Coupon4

Probabilistic model of Coupon Collector’s program for four coupons, taken from [18].

⬇

f := 0

g := 0

a := 0

b := 0

c := 0

d := 0

while true:

f := 1 [1/2] 0

g := 1 [1/2] 0

a := a + (1-a)fg

b := b + (1-b)f(1-g)

c := c + (1-c)*(1-f)*g

d := d + (1-d)(1-f)(1-g)

5.3 Random_walk_1D_cts

A variation of random walk in one dimension with assignments from continuous distributions taken from [18].

⬇

v := 0

x := 0

while true:

v := u(0,1)

x := x + v [7/10] x - v

5.4 Sum_rnd_series

A program modeling Sum of Random Series game taken from [6].

⬇

n := 0

x := 0

while true:

n := n + 1

x := x + n [1/2] x

5.5 Product_dep_var

A program modeling Product of Dependent Random Variables game taken from [6].

⬇

f := 0

x := 0

y := 0

p := 0

while true:

f := 0 [1/2] 1

x := x + f

y := y + 1 - f

p := x*y

5.6 Random_walk_2D

A variation of random walk in two dimension as in [4, 18]. Each direction is chosen with equal probability.

⬇

h := 0

x := 0

y := 0

while true:

h := 1 [1/2] 0

x := x-h [1/2] x +h

y := y+(1-h) [1/2] y-(1-h)

5.7 Binomial

Another classic example, modeling binomial distribution. Appeared also in [6, 8, 14]. However, we consider the program to be parametric, computing invariants for arbitrary values of $p$.

⬇

x := 0

while true:

x := x + 1 [p] x

5.8 StutteringA

Program 1(A) from Introduction.

5.9 StutteringB

Program 1(B) from Introduction.

5.10 StutteringC

Program 1(C) from Introduction.

5.11 StutteringD

Program 1(D) from Introduction.

5.12 StutteringP

A variation of program 1(A) from Introduction with $d=1$, parametrized w.r.t. $p$.

⬇

f := 0

x := -1

y := 1

s := 0

while true:

f := 1 [p] 0

x := x + f*u(0,2)

y := y + f*u(0,4)

s := x + y

5.13 Square

Our own program with polynomial assignments.

⬇

x := 0; y := 1

while true:

x := x+2 [1/2] x

y := x^2

Experimental Results with Moment-Based Invariants.

Results of our evaluation are presented in Table 1. While Algorithm 1 can compute invariants over arbitrary $k$th higher-order moments, due to lack of space and readability, Table 1 lists only our moment-based invariants up to the third moment (i.e. $k\leq 3$), that is for expected values, second- and third-order moments. The first column of Table 1 lists the benchmark name, whereas the second column gives the degree of the moments (i.e. $k=1,2,3$) for which we compute invariants. The third column reports the timings (in seconds) our implementation needed to derive invariants. The last column shows our moment-based invariants; for readability, we decided to omit intermediary invariants (up to $30$ for some programs) and only show the most relevant invariants.

We could not perform a fair practical comparison with other existing methods: to the best of our knowledge, existing works, such as [14, 11, 2, 18], require user guidance/templates/hints. Further, existing techniques do not support symbolic probabilities and/or parametrised distributions - which are, for example, required in the analysis of programs StutteringA, StutteringC, StutteringP of Table 1. We also note that examples Coupon, StutteringC, StutteringP involve non-linear probabilistic updates hindering automation in existing methods, while such updates can naturally be encoded as moment-based recurrences in our framework. We finally note that while second-order moments are computed only by [18], but with the help of user-provided templates, no existing approaches compute moments for $k\geq 3$. Our experiments show that inferring third-order moments are in general not expensive; yet, for examples StutteringA, StutteringC, StutteringP with parametrized distributions/probabilities more computation time is needed.

6 Related Work

Despite the impressive recent advancements, probabilistic model checking tools [1] (e.g., PRISM [19], STORM [7] and MRMC [15]) are not able to handle programs with unbounded and real variables. Model checking algorithms suffer from the state explosion problem and their performance in terms of time and memory consumption degrades as the number of reachable states to be considered increases. Furthermore, probabilistic model checking tools have no support for invariant generation. Our approach, based on symbolic summation over probabilistic expressions, can instead analyse probabilistic programs with a potentially infinite number of reachable states.

In [21], one of the first deductive frameworks to reason about probabilistic programs was proposed by annotating probabilistic programs with real-valued expressions over the expected values of program variables. Of particular interest are the annotations as quantitative invariants, summarising loop behaviors. The setting of [21] considers probabilistic programs where the stochastic inputs are restricted to discrete distributions with finite support and can deal also with demonic non-deterministic choice. Although our approach does not yet support demonic non-determinism, we are not restricted to discrete input distributions as long as we know their moments (e.g., the Gaussian distribution is characterised only by two moments: the mean and the variance). Moreover, our work is not restricted to quantitative invariants as invariants over expected values of program variables. Rather, we generate moment-based invariants that precisely capture invariant properties of higher-order and mixed moments of program variables.

Katoen et al. provided in [14] the first semi-automatic and complete method synthesising the linear quantitative invariants needed by [21]. The work of [14], implemented in PRINSYS [11], consists in annotating a loop with a linear template invariants and uses a constraint solver to find the parameters for which the template yields an invariant. The works [8, 6] also synthesize non-linear quantitative invariants.

Another related line of research is given in [2], where martingales are used to compute invariants of probabilistic programs. The martingales generated by [2] however heavily depend on the user-provided hints and hence less generic hints yield less expressive/precise invariants. Moreover, of [2] mainly focuses on invariants over expected values and it remains unclear which extensions of martingales need to be considered to compute higher-order moments. The work of [18] addresses such generalizations of martingales for computing higher-order moments of program variables, with the overall goal of approximating runtimes of randomized programs. The approach in [18] is however again restricted to user-provided templates. Unlike the works of [14, 11, 8, 6, 2, 18], our work does not rely on a priori given templates/hints, but computes the most precise invariant expression over higher-order or mixed moments of program variables. To do so, we use symbolic summation to compute closed forms of higher-order moments. In addition, Prob-solvable loops support parametrized distributions and symbolic probabilities, which is not the case of [2, 18].

There are two orthogonal problems related to quantitative invariants generation: program termination [22, 9] and worst-case execution [3, 13, 5]. The first is to assess whether a probabilistic program terminates with probability 1 or if the expected time of termination is bounded. In principle, one can use our approach to solve this class of problems for Prob-solvable loops, but this is not the focus of this paper. The second class of problems is related to finding bounds over the expected values. In [3] the authors consider bounds also over higher-order moments for a specific class of probabilistic programs with probabilistic affine assignments. This approach can handle also nonlinear terms using interval arithmetics and fresh variables, at the price to produce very conservative bounds. On the contrary our approach supports natively probabilistic polynomial assignments (in the form of Prob-solvable loops) and provides a precise symbolic expression over higher-order moments.

7 Conclusion

We introduced a novel approach for automatically generating moment-based invariants of a subclass of probabilistic programs (PPs), called Prob-solvable loops, with polynomial assignments over random variables and parametrised distributions. We combine methods from symbolic summation and statistics to derive invariants over higher-order moments, such as expected values or variances, of program variables. To the best of our knowledge, our approach is the first method computing higher-order moments of PPs fully automatically and the first to handle PPs with parametrised distributions.