Central Moment Analysis for Cost Accumulators in Probabilistic Programs
Technical Report

Di Wang

Carnegie Mellon University

Jan Hoffmann

Carnegie Mellon University

Thomas Reps

University of Wisconsin

For probabilistic programs, it is usually not possible to automatically
derive exact information about their properties, such as the distribution
of states at a given program point.
Instead, one can attempt to derive approximations, such as upper bounds
on tail probabilities.

Such bounds can be obtained via concentration inequalities, which rely
on the moments of a distribution, such as
the expectation (the first raw moment) or
the variance (the second central moment).

Tail bounds obtained using central moments are often
tighter than the ones obtained using raw moments, but
automatically analyzing central moments is more challenging.

This paper presents an analysis for probabilistic programs that
automatically derives symbolic upper and lower bounds on
variances, as well as higher central moments, of cost accumulators.

To overcome the challenges of higher-moment analysis, it generalizes
analyses for expectations with an algebraic abstraction that
simultaneously analyzes different moments, utilizing relations between

A key innovation is the notion of moment-polymorphic recursion, and a
practical derivation system that handles recursive functions.

The analysis has been implemented using a template-based technique
that reduces the inference of polynomial bounds to linear programming.

Experiments with our prototype central-moment analyzer show that,
despite the analyzer's upper/lower bounds on various quantities,
it obtains tighter tail bounds than an existing system that uses only raw
moments, such as expectations.

<concept_desc>Theory of computation Probabilistic computation</concept_desc>
<concept_desc>Theory of computation Program analysis</concept_desc>
<concept_desc>Theory of computation Automated reasoning</concept_desc>

[500]Theory of computation Probabilistic computation
[500]Theory of computation Program analysis
[500]Theory of computation Automated reasoning

§ INTRODUCTION

Probabilistic programs <cit.> can be used to manipulate uncertain quantities modeled by probability distributions and random control flows.

Uncertainty arises naturally in Bayesian networks, which capture statistical dependencies (e.g., for diagnosis of diseases <cit.>), and cyber-physical systems, which are subject to sensor errors and peripheral disturbances (e.g., airborne collision-avoidance systems <cit.>).

Researchers have used probabilistic programs to implement and analyze
randomized algorithms <cit.>, cryptographic protocols <cit.>, and machine-learning algorithms <cit.>.

A probabilistic program propagates uncertainty through a computation, and produces a distribution over results.

In general, it is not tractable to compute the result distributions of probabilistic programs automatically and precisely: Composing simple distributions can quickly complicate the result distribution, and randomness in the control flow can easily lead to state-space explosion.

Monte-Carlo simulation <cit.> is a common approach to study the result distributions, but the technique does not provide formal guarantees, and can sometimes be inefficient <cit.>.

In this paper, we focus on a specific yet important kind of uncertain quantity: cost accumulators, which are program variables that can only be incremented or decremented through the program execution.

Examples of cost accumulators include
termination time <cit.>,
rewards in Markov decision processes (MDPs) <cit.>,
position information in control systems <cit.>,
and cash flow during bitcoin mining <cit.>.

Recent work <cit.> has proposed successful static-analysis approaches that leverage
aggregate information of a cost accumulator $X$, such as $X$'s expected value $\expe[X]$ (i.e., $X$'s “first moment”).

The intuition why it is beneficial to compute aggregate information—in lieu of distributions—is that aggregate measures like expectations
abstract distributions to a single number, while still indicating non-trivial properties.

Moreover, expectations are transformed by statements in a probabilistic program in a manner similar to the weakest-precondition transformation of formulas in a non-probabilistic program <cit.>.

One important kind of aggregate information is moments.

In this paper, we show how to obtain central moments (i.e., $\expe[(X-\expe[X])^k]$ for any $k \ge 2$),
whereas most previous work focused on raw moments (i.e., $\expe[X^k]$ for any $k \ge 1$).

Central moments can provide more information about distributions.

For example, the variance $\vari[X]$ (i.e., $\expe[(X-\expe[X])^2]$, $X$'s “second central moment”) indicates how $X$ can deviate from its mean,
the skewness (i.e., $\frac{ \expe[(X-\expe[X])^3] }{ (\vari[X])^{\sfrac{3}{2}} }$, $X$'s “third standardized moment”) indicates how lopsided the distribution of $X$ is,
and the kurtosis (i.e., $\frac{ \expe[(X - \expe[X])^4 ] }{ (\vari[X])^2 }$, $X$'s “fourth standardized moment”) measures the heaviness of the tails of the distribution of $X$.

One application of moments is to answer queries about tail bounds, e.g., the assertions about probabilities of the form $\prob[X \ge d]$, via concentration-of-measure inequalities from probability theory <cit.>.

With central moments, we find an opportunity to obtain more precise tail bounds of the form $\prob[X \ge d]$, and become able to derive bounds on tail probabilities of the form $\prob[|X - \expe[X]| \ge d]$.

Central moments $\expe[(X-\expe[X])^k]$ can be seen as polynomials of raw moments $\expe[X],\cdots,\expe[X^k]$, e.g., the variance $\vari[X] = \expe[(X-\expe[X])^2]$ can be rewritten as $\expe[X^2]-\expe^2[X]$, where $\expe^k[X]$ denotes $(\expe[X])^k$.

To derive bounds on central moments, we need both upper and lower bounds on the raw moments, because of the presence of subtraction.

For example, to upper-bound $\vari[X]$, a static analyzer needs to have an upper bound on $\expe[X^2]$ and a lower bound on $\expe^2[X]$.

In this work, we present and implement the first fully automatic analysis for deriving symbolic interval
bounds on higher central moments for cost accumulators in probabilistic programs with general recursion
and continuous distributions.

One challenge is to support interprocedural reasoning to reuse analysis results for functions.

Our solution makes use of a “lifting” technique from the natural-language-processing community.

That technique derives an algebra for second moments from an algebra for first moments <cit.>.

We generalize the technique to develop moment semirings,
and use them to derive a novel frame rule to handle function calls with moment-polymorphic recursion (see <ref>).

Recent work has successfully automated inference of upper <cit.> or lower bounds <cit.> on the expected cost of probabilistic programs.

<cit.> developed a system to derive upper bounds on higher raw moments of program runtimes.

However, even in combination, existing approaches cannot solve tasks such as
deriving a lower bound on the second raw moment of runtimes, or
deriving an upper bound on the variance of accumulators that count live heap cells.

<ref>(a) summarizes the features of related work on moment inference for probabilistic programs.

To the best of our knowledge, our work is the first moment-analysis tool that supports all of the listed programming and analysis features.

<ref>(b) and (c) compare our work with related work in terms of tail-bound analysis
on a concrete program (see <ref>).

The bounds are derived for the cost accumulator 𝕀tick in a random-walk program that we will present in <ref>.

It can be observed that for $d \ge 20$, the most precise tail bound for 𝕀tick is the one obtained via an upper bound on the variance $\vari[\id{tick}]$ ($\id{tick}$'s second central moment).

feature     <cit.>     <cit.>     <cit.>     <cit.>     this work

continuous distributions

non-monotone costs

higher moments

interval bounds

<cit.>      <cit.>     this work
Derived      2*$\expe[\id{tick}] \le 2d + 4$     $\expe[\id{tick}^2] \le {}$     $\vari[\id{tick}] \le {} $

bound         $4d^2+22d+28$     $22d+28$
Moment type     raw     raw     central
Concentration     Markov     Markov     2*Cantelli

inequality     ($\text{degree} = 1$)     ($\text{degree} = 2$)
Tail bound      2* $\approx \frac{1}{2}$     2*$\approx \frac{1}{4} $     2*$\xrightarrow{d \to \infty} 0 $

$\prob[\id{tick} \ge 4d]$

(a) Comparison in terms of supporting features. (b) Comparison in terms of moment bounds for the running example. (c) Comparison in terms of derived tail bounds.

Our work incorporates ideas known from the literature:

* Using the expected-potential method (or ranking super-martingales) to derive upper bounds on the expected program runtimes or monotone costs <cit.>.

* Using the Optional Stopping Theorem from probability theory to ensure the soundness of lower-bound inference for probabilistic programs <cit.>.

* Using linear programming (LP) to efficiently automate the (expected) potential method for (expected) cost analysis <cit.>.

The contributions of our work are as follows:

* We develop moment semirings to compose the moments for a cost accumulator from two computations, and to enable interprocedural reasoning about higher moments.

* We instantiate moment semirings with the symbolic interval domain, use that to develop a derivation system for interval bounds on higher central moments for cost accumulators, and automate the derivation via LP solving.

* We prove the soundness of our derivation system for programs that satisfy the criterion of our recent extension to the Optional Stopping Theorem,
and develop an algorithm for checking this criterion automatically.

We implemented our analysis
and evaluated it

on a broad suite of benchmarks from the literature.

Our experimental results show that on a variety of examples, our analyzer is able to use higher central moments to obtain tighter tail bounds on program runtimes than the system of <cit.>, which uses only upper bounds on raw moments.

§ OVERVIEW

In this section, we demonstrate the expected-potential method for both first-moment analysis (previous work) and higher central-moment analysis (this work) (<ref>),

and discuss the challenges to supporting interprocedural reasoning and to ensuring the soundness of our approach (<ref>).

The program in <ref> implements a bounded, biased random walk.

The main function consists of a single statement “$\iinvoke{\mathsf{rdwalk}}$” that invokes a recursive function.

The variables $x$ and $d$ represent the current position and the ending position of the random walk, respectively.

We assume that $d>0$ holds initially.

In each step, the program samples the length of the current move from a uniform distribution on the interval $[{-1},2]$.

The statement $\itick{1}$ adds one to a cost accumulator that counts the number of steps before
the random walk ends.

We denote this accumulator by 𝕀tick in the rest of this section.

The program terminates with probability one and its expected accumulated cost is bounded by

func $\mathsf{rdwalk}$() begin
${ \color{ACMDarkBlue} \{ \; 2(d-x)+4 \;\} }$

if $\cbin{<}{x}{d}$ then
${\color{ACMDarkBlue} \{\;  2(d-x)+4 \;\}  }$

$\isample{t}{ \kw{uniform}({-1},2)}$;

${\color{ACMDarkBlue} \{\;  2(d - x - t) + 5  \;\} }$

$ \iassign{x}{  \eadd{x}{t} }$;

${\color{ACMDarkBlue} \{\;  2(d - x) + 5 \;\} }$

call $\mathsf{rdwalk}$;

${\color{ACMDarkBlue} \{ \; 1 \; \} }$

${\color{ACMDarkBlue} \{ \; 0 \; \} } $

# $\mathsf{XID} \defeq \{x,d,t\}$

# $\mathsf{FID} \defeq \{\mathsf{rdwalk} \}$

# pre-condition: $\{d > 0 \}$

func $\mathsf{main}$() begin

call $\mathsf{rdwalk}$

A bounded, biased random walk, implemented using recursion. The annotations show the derivation of an upper bound on the expected accumulated cost.

§.§ The Expected-Potential Method for Higher-Moment Analysis

Our approach to higher-moment analysis is inspired by the expected-potential method <cit.>, which is also known as ranking super-martingales <cit.>, for expected-cost bound analysis of probabilistic programs.

The classic potential method of amortized analysis <cit.> can be automated to derive symbolic cost bounds for non-probabilistic programs <cit.>.

The basic idea is to define a potential function $\phi: \Sigma \to \bbR^+$ that maps program states $\sigma \in \Sigma$ to nonnegative numbers,
where we assume each state $\sigma$ contains a cost-accumulator component $\sigma.\alpha$.

If a program executes with initial state $\sigma$ to final state $\sigma'$,
then it holds that $\phi(\sigma) \ge (\sigma'.\alpha - \sigma.\alpha) + \phi(\sigma')$,
where $(\sigma'.\alpha - \sigma.\alpha)$ describes the accumulated cost from $\sigma$ to $\sigma'$.

The potential method also enables compositional reasoning: if a statement $S_1$ executes from $\sigma$ to $\sigma'$ and a statement $S_2$ executes from $\sigma'$ to $\sigma''$, then we have $\phi(\sigma) \ge (\sigma'.\alpha - \sigma.\alpha) + \phi(\sigma')$ and $\phi(\sigma') \ge (\sigma''.\alpha - \sigma'.\alpha) + \phi(\sigma'')$; therefore, we derive $\phi(\sigma) \ge (\sigma''.\alpha - \sigma.\alpha) + \phi(\sigma'')$ for the sequential composition $S_1;S_2$.

For non-probabilistic programs, the initial potential provides an upper bound on the accumulated cost.

This approach has been adapted to reason about expected costs of probabilistic programs <cit.>.

To derive upper bounds on the expected accumulated cost of a program $S$ with initial state $\sigma$, one needs to take into consideration the distribution of all possible executions.

More precisely, the potential function should satisfy the following property:

\begin{equation}\label{Eq:ExpectedPotentialInequality}
\phi(\sigma) \ge \expe_{\sigma' \sim \interp{S}(\sigma)}[ C(\sigma,\sigma') + \phi(\sigma')  ],
\end{equation}

where the notation $\expe_{x \sim \mu}[f(x)]$ represents the expected value of $f(x)$, where $x$ is drawn from the distribution $\mu$, $\interp{S}(\sigma)$ is the distribution over final states of executing $S$ from $\sigma$, and $C(\sigma,\sigma') \defeq \sigma'.\alpha - \sigma.\alpha$ is the execution cost from $\sigma$ to $\sigma'$.

<ref> annotates the $\mathsf{rdwalk}$ function from <ref> with the derivation of an upper bound on the expected accumulated cost.

The annotations, taken together, define an expected-potential function $\phi : \Sigma \to \bbR^+$ where a program state $\sigma \in \Sigma$ consists of a program point and a valuation for program variables.

To justify the upper bound $2(d-x)+4$ for the function $\mathsf{rdwalk}$, one has to show that the potential right before the $\itick{1}$ statement should be at least $1$.

This property is established by backward reasoning on the function body:

* For $\iinvoke{\mathsf{rdwalk}}$, we apply the “induction hypothesis” that the expected cost of the function $\mathsf{rdwalk}$ can be upper-bounded by $2(d-x)+4$.

Adding the $1$ unit of potential need by the tick statement, we obtain $2(d-x)+5$ as the pre-annotation of the function call.

* For $\iassign{x}{\eadd{x}{t}}$, we substitute $x$ with $\eadd{x}{t}$ in the post-annotation of this statement to obtain the pre-annotation.

* For $\isample{t}{\kw{uniform}({-1},2)}$, because its post-annotation is $2(d-x-t)+5$, we compute its pre-annotation as

\begin{aligned}
& \expe_{t \sim \kw{uniform}({-1},2)}[2(d-x-t)+5] \\
={} & 2(d-x)+5 - 2 \cdot \expe_{t \sim \kw{uniform}({-1},2)}[t] \\
={} & 2(d-x)+5-2 \cdot \sfrac{1}{2} = 2(d-x)+4,
\end{aligned}

which is exactly the upper bound we want to justify.

Our approach

In this paper, we focus on derivation of higher central moments.

Observing that a central moment $\expe[(X-\expe[X])^k]$ can be rewritten as a polynomial of raw moments $\expe[X],\cdots,\expe[X^k]$,
we reduce the problem of bounding central moments to reasoning about upper and lower bounds on raw moments.

For example, the variance can be written as $\vari[X] = \expe[X^2] - \expe^2[X]$, so it suffices to analyze the upper bound of
the second moment $\expe[X^2]$ and the lower bound on the square of the first moment $\expe^2[X]$.

For higher central moments, this approach requires both upper and lower bounds on higher raw moments.

For example, consider the fourth central moment of a nonnegative random variable $X$:

\expe[(X-\expe[X])^4] = \expe[X^4]-4\expe[X^3]\expe[X]+6\expe[X^2]\expe^2[X] - 3\expe^4[X].

Deriving an upper bound on the fourth central moment requires lower bounds on the first ($\expe[X]$) and third ($\expe[X^3]$) raw moments.

We now sketch the development of moment semirings.

We first consider only the upper bounds on higher moments of nonnegative costs.

To do so, we extend the range of the expected-potential function $\phi$ to real-valued vectors $(\bbR^+)^{m+1}$, where $m \in \bbN$ is the degree of the target moment.

We update the potential inequality (<ref>) as follows:

\begin{equation}\label{Eq:OverApproxIntervalPotentialInequality}
\phi(\sigma) \ge \expe_{\sigma' \sim \interp{S}(\sigma)} [ \many{\langle C(\sigma,\sigma')^k\rangle_{0\le k\le m}} \otimes \phi(\sigma')    ],
\end{equation}

where $\many{\tuple{v_k}_{0 \le k \le m}}$ denotes an $(m+1)$-dimensional vector, the order $\le$ on vectors is defined pointwise, and $\otimes$ is some composition operator.

Recall that $\interp{S}(\sigma)$ denotes the distribution over final states of executing $S$ from $\sigma$, and $C(\sigma,\sigma')$ describes the cost for the execution from $\sigma$ to $\sigma'$.

Intuitively, for $\phi(\sigma) = \many{ \tuple{\phi(\sigma)_k}_{0 \le k  \le m}}$ and each $k$, the component $\phi(\sigma)_k$ is an upper bound on the $k$-th moment of the cost for the computation starting from $\sigma$.

The $0$-th moment is the termination probability of the computation, and we assume it is always one for now.

We cannot simply define $\otimes$ as pointwise addition because, for example, $(a+b)^2 \neq a^2+b^2$ in general.

If we think of $b$ as the cost for some probabilistic computation, and we prepend a constant cost $a$ to the computation, then by linearity of expectations, we have $\expe[(a+b)^2] = \expe[a^2+2ab+ b^2] = a^2 + 2 \cdot a \cdot \expe[b] + \expe[b^2]$, i.e., reasoning about the second moment requires us to keep track of the first moment.

Similarly, we should have

\begin{equation*}
\phi(\sigma)_2 \ge \expe_{\sigma' \sim \interp{S}(\sigma)}[ C(\sigma,\sigma')^2 + 2 \cdot C(\sigma,\sigma') \cdot \phi(\sigma')_1 + \phi(\sigma')_2 ],
\end{equation*}

for the second-moment component,
where $\phi(\sigma')_1$ and $\phi(\sigma')_2$ denote $\expe[b]$ and $\expe[b^2]$, respectively.

Therefore, the composition operator $\otimes$ for second-moment analysis (i.e., $m=2$) should be defined as

\begin{equation}\label{Eq:SecondMomentComposition}
\tuple{ 1, r_1,s_1 } \otimes \tuple{ 1, r_2,s_2 } \defeq \tuple{ 1, r_1 + r_2, s_1 + 2r_1r_2 + s_2 }.
\end{equation}

func $\mathsf{rdwalk}$() begin
${\color{ACMDarkBlue} \{ \; \langle 1, \; 2(d-x) + 4, \; 4(d - x)^2 + 22(d - x) + 28 \rangle \; \} }$

if $\cbin{<}{x}{d}$ then
${\color{ACMDarkBlue} \{\; \langle 1, \; 2(d - x) + 4, \; 4(d - x)^2 + 22(d -x ) + 28 \rangle \; \}  }$

$\isample{t}{ \kw{uniform}({-1},2)}$;

${\color{ACMDarkBlue} \{  \; \langle 1, \; 2(d-x-t)+5, \; 4(d-x-t)^2 + 26(d - x - t) + 37 \rangle \;  \} }$

$ \iassign{x}{  \eadd{x}{t} }$;

${\color{ACMDarkBlue} \{\; \langle 1, \; 2(d - x) + 5, \; 4(d  - x)^2 + 26(d-x) + 37 \rangle \; \} }$

call $\mathsf{rdwalk}$;

${\color{ACMDarkBlue} \{  \; \tuple{1, \; 1, \; 1}  \; \} }$

${\color{ACMDarkBlue} \{ \; \tuple{1, \; 0, \; 0} \; \} } $

Derivation of an upper bound on the first and second moment of the accumulated cost.

<ref> annotates the $\mathsf{rdwalk}$ function from <ref> with the derivation of an upper bound on both the first and second moment of the accumulated cost.

To justify the first and second moment of the accumulated cost for the function $\mathsf{rdwalk}$,

we again perform backward reasoning:

* For $\itick{1}$, it transforms a post-annotation $a$ by $\lambda a. (\tuple{1,1,1} \!\otimes\! a)$; thus, the pre-annotation is $\tuple{1,1,1} \otimes \tuple{1,0,0} = \tuple{1,1,1}$.

* For $\iinvoke{\mathsf{rdwalk}}$, we apply the “induction hypothesis”, i.e., the upper bound shown on line 2. We use the $\otimes$ operator to compose the induction hypothesis with the post-annotation of this function call:

\begin{aligned}
& \tuple{ 1, 2(d\!-\!x)\!+\!4 ,4(d\!-\!x)^2\!+\!22(d\!-\!x)\!+\!28 } \otimes \tuple{1,1,1} \\[-3pt]
={} & \tuple{1, 2(d\!-\!x)\!+\!5, (4(d\!-\!x)^2\!+\!22(d\!-\!x)\!+\!28)\!+\!2\!\cdot\! (2(d\!-\!x)\!+\!4) \!+\! 1  } \\[-3pt]
={} & \tuple{1, 2(d\!-\!x)\!+\!5, 4(d\!-\!x)^2\!+\!26(d\!-\!x)\!+\!37 }.
\end{aligned}

* For $\iassign{x}{\eadd{x}{t}}$, we substitute $x$ with $x+t$ in the post-annotation of this statement to obtain the pre-annotation.

* For $\isample{t}{\kw{uniform}({-1},2)}$, because the post-annotation involves both $t$ and $t^2$, we compute from the definition of uniform distributions that

\[
\expe_{t \sim \kw{uniform}({-1},2)}[t]=\sfrac{1}{2}, \quad \expe_{t \sim \kw{uniform}({-1},2)}[t^2] = 1.
\]

Then the upper bound on the second moment is derived as follows:

\begin{aligned}
& \expe_{t \sim \kw{uniform}({-1},2)}[4(d\!-\!x\!-\!t)^2\!+\!26(d\!-\!x\!-\!t)\!+\!37] \\
={} & (4(d\!-\!x)^2\!+\!26(d\!-\!x)\!+\!37)  \!-\! (8(d\!-\!x)\!+\!26) \!\cdot\! \expe_{t \sim \kw{uniform}({-1},2)}[t ] \\[-3pt]
& + 4 \!\cdot\! \expe_{t \sim \kw{uniform}({-1},2)}[t^2] \\
={} & 4(d\!-\!x)^2 \!+\!22(d\!-\!x) \!+\!28,
\end{aligned}

which is the same as the desired upper bound on the second moment of the accumulated cost for the function $\mathsf{rdwalk}$. (See <ref>, line 2.)

We generalize the composition operator $\otimes$ to moments with arbitrarily high degrees, via a family of algebraic structures, which we name moment semirings (see <ref>).

These semirings are algebraic in the sense that they can be instantiated with any partially ordered semiring, not just $\bbR^+$.

Interval bounds

Moment semirings not only provide a general method to analyze higher moments, but also enable reasoning about upper and lower bounds on moments simultaneously.

The simultaneous treatment is also essential for analyzing programs with non-monotone costs (see <ref>).

We instantiate moment semirings with the standard interval semiring $\calI = \{ [a,b] \mid a \le b \}$.

The algebraic approach allows us to systematically incorporate the interval-valued bounds, by reinterpreting operations in <ref> under $\calI$:

\begin{aligned}
& \tuple{ [1,1], [r_1^{\m{L}},r_1^{\m{U}}], [s_1^{\m{L}},s_1^{\m{U}}] } \otimes \tuple{ [1,1] , [r_2^{\m{L}},r_2^{\m{U}}], [s_2^{\m{L}},s_2^{\m{U}}] } \\
{} \defeq {} & \langle [1,1], [r_1^{\m{L}},r_1^{\m{U}}] +_\calI  [r_2^{\m{L}},r_2^{\m{U}}], \\[-3pt]
& \hphantom{\langle} [s_1^{\m{L}},s_2^{\m{U}}]+_\calI 2 \cdot ([r_1^{\m{L}},r_1^{\m{U}}] \cdot_\calI [r_2^{\m{L}},r_2^{\m{U}}]) +_\calI [s_2^{\m{L}},s_2^{\m{U}}] \rangle \\
={} & \langle [1,1], [r_1^{\m{L}}+r_2^{\m{L}},r_1^{\m{U}}+r_2^{\m{U}}], [s_1^{\m{L}}+2 \cdot \min S+s_2^{\m{L}}, s_1^{\m{U}}+2 \cdot \max S+s_2^{\m{U}}] \rangle ,
\end{aligned}

where $S \defeq \{ r_1^{\m{L}}r_2^{\m{L}}, r_1^{\m{L}}r_2^{\m{U}}, r_1^{\m{U}}r_2^{\m{L}}, r_1^{\m{U}}r_2^{\m{U}} \}$.

We then update the potential inequality <ref> as follows:

\begin{equation*}%\label{Eq:IntervalApproxIntervalPotentialInequality}
\phi(\sigma) \sqsupseteq \expe_{\sigma' \sim \interp{S}(\sigma)} [ \many{\langle [C(\sigma,\sigma')^k,C(\sigma,\sigma')^k] \rangle_{0\le k \le m}} \otimes \phi(\sigma')    ],
\end{equation*}

where the order $\sqsubseteq$ is defined as pointwise interval inclusion.

Suppose that the interval bound on the first moment of the accumulated cost of the $\m{rdwalk}$ function from <ref> is $[2(d-x),2(d-x)+4]$.

We can now derive the upper bound on the variance $\vari[\id{tick}] \le 22d+28$ shown in <ref>(b) (where we substitute $x$ with $0$ because the main function initializes $x$ to $0$ on line 5 in <ref>):

\begin{aligned}
\vari[\id{tick}] & = \expe[\id{tick}^2] - \expe^2[\id{tick}] \\[-3pt]
& \le (\text{upper bnd. on $\expe[\id{tick}^2]$}) \!- \! (\text{lower bnd. on $\expe[\id{tick}]$})^2 \\[-3pt]
& = (4d^2+22d+28) - (2d)^2 = 22d+28.
\end{aligned}

In <ref>, we describe how we use moment bounds to derive the tail bounds shown in <ref>(c).

§.§ Two Major Challenges

Interprocedural reasoning

Recall that in the derivation of <ref>, we use the $\otimes$ operator to compose the upper bounds on moments for $\iinvoke{\m{rdwalk}}$ and its post-annotation $\tuple{1,1,1}$.

However, this approach does not work in general, because the post-annotation might be symbolic (e.g., $\tuple{1,x,x^2}$) and the callee might mutate referenced program variables (e.g., $x$).

One workaround is to derive a pre-annotation for each possible post-annotation of a recursive function, i.e., the moment annotations for a recursive function is polymorphic.

This workaround would not be effective for non-tail-recursive functions:
for example, we need to reason about the $\m{rdwalk}$ function in <ref> with infinitely many post-annotations $\tuple{1,0,0}$, $\tuple{1,1,1}$, $\tuple{1,2,4}$, …, i.e., $\tuple{1,i,i^2}$ for all $i \in \bbZ^+$.

Our solution to moment-polymorphic recursion is to
introduce a combination operator $\oplus$ in a way that
if $\phi_1$ and $\phi_2$ are two expected-potential functions, then

\[
\phi_1(\sigma) \oplus \phi_2(\sigma) \!\ge\! \expe_{\sigma'  \!\sim\! \interp{S}(\sigma)}[ \many{ \tuple{C(\sigma,\sigma')^k}_{0 \!\le\! k \!\le\! m} } \otimes (\phi_1(\sigma') \oplus \phi_2(\sigma')) ].
\]

We then use the $\oplus$ operator to derive a frame rule:
{  Q_1   }    S    {   Q_1'   }

{  Q_2   }    S    {   Q_2'   }

{  Q_1 ⊕Q_2   }    S    {   Q_1' ⊕Q_2'   }

We define $\oplus$ as pointwise addition, i.e., for second moments,

\begin{equation}\label{Eq:SndMomentCombine}
\tuple{p_1,r_1,s_1} \oplus \tuple{p_2,r_2,s_2} \defeq \tuple{p_1+p_2,r_1+r_2,s_1+s_2},
\end{equation}

and because the $0$-th-moment (i.e., termination-probability) component is no longer guaranteed to be one, we redefine $\otimes$ to consider the termination probabilities:

\begin{equation}\label{Eq:SndMomentExtend}
\tuple{p_1,r_1,s_1} \otimes \tuple{p_2,r_2,s_2} \defeq \tuple{p_1 p_2, p_2r_1+p_1r_2, p_2s_1+2r_1r_2 + p_1s_2}.
\end{equation}

As we will show in <ref>,
the composition operator $\otimes$ and combination operator $\oplus$ form a moment semiring;
consequently, we can use algebraic properties of semirings (e.g., distributivity) to aid higher-moment analysis.

For example, a vector $\tuple{0,r_1,s_1}$ whose termination-probability component is zero does not seem to make sense, because moments with respect to a zero distribution should also be zero.

However, by distributivity, we have

\begin{aligned}
& \tuple{1,r_3,s_3} \otimes \tuple{1,r_1+r_2,s_1+s_2} \\
={} & \tuple{1,r_3,s_3} \otimes (\tuple{0,r_1,s_1} \oplus \tuple{1,r_2,s_2}) \\
={} & (\tuple{1,r_3,s_3} \otimes \tuple{0,r_1,s_1}) \oplus ( \tuple{1,r_3,s_3} \oplus \tuple{1,r_2,s_2}).
\end{aligned}

If we think of $\tuple{1,r_1+r_2,s_1+s_2}$ as a post-annotation of a computation whose moments are bounded by $\tuple{1,r_3,s_3}$,
the equation above indicates that we can use $\oplus$ to decompose the post-annotation into subparts, and then reason about each subpart separately.

This fact inspires us to develop a decomposition technique for moment-polymorphic recursion.

With the $\oplus$ operator and the frame rule, we only need to analyze the $\m{rdwalk}$ function from <ref> with three post-annotations: $\tuple{1,0,0}$, $\tuple{0,1,1}$, and $\tuple{0,0,2}$, which form a kind of “elimination sequence.”

We construct this sequence in an on-demand manner;
the first post-annotation is the identity element $\tuple{1,0,0}$ of the moment semiring.

For post-annotation $\tuple{1,0,0}$, as shown in <ref>, we need to know the moment bound for $\m{rdwalk}$ with the post-annotation $\tuple{1,1,1}$.

Instead of reanalyzing $\m{rdwalk}$ with the post-annotation $\tuple{1,1,1}$,
we use the $\oplus$ operator to compute the “difference” between it and the previous post-annotation $\tuple{1,0,0}$.

Observing that $\tuple{1,1,1} = \tuple{1,0,0} \oplus \tuple{0,1,1}$, we now analyze $\m{rdwalk}$ with $\tuple{0,1,1}$ as the post-annotation:

$\iinvoke{\m{rdwalk}}$;  $\{ \;  \tuple{0, \; 1, \; 3} \;  \}$    # $=\tuple{1, \; 1, \; 1} \otimes \tuple{0, \; 1, \; 1} $

$\itick{1}$  $\{ \; \tuple{0,\; 1, \;1} \; \}$

Again, because $\tuple{0,1,3} = \tuple{0,1,1} \oplus \tuple{0,0,2}$, we need to further analyze $\m{rdwalk}$ with $\tuple{0,0,2}$ as the post-annotation:

$\iinvoke{\m{rdwalk}}$;  $\{ \;  \tuple{0, \; 0, \; 2} \;  \}$    # $=\tuple{1, \; 1, \; 1} \otimes \tuple{0, \; 0, \; 2}$

$\itick{1}$  $\{ \; \tuple{0,\; 0, \;2} \; \}$

With the post-annotation $\tuple{0,0,2}$, we can now reason monomorphically without analyzing any new post-annotation!

We can perform a succession of reasoning steps similar to what we have done in <ref> to justify the following bounds (“unwinding” the elimination sequence):

* ${\color{ACMDarkBlue} \{  \tuple{ 0,  0 , 2 }  \} } ~ {\m{rdwalk}} ~ { \color{ACMDarkBlue} \{  \tuple{0,  0,  2}  \} }$:
Directly by backward reasoning with the post-annotation $\tuple{0,0,2}$.

* ${\color{ACMDarkBlue} \{  \tuple{ 0,  1 , 4(d-x) + 9  }  \} } ~ {\m{rdwalk}} ~ { \color{ACMDarkBlue} \{  \tuple{0,  1, 1}  \} }$:
To analyze the recursive call with post-annotation $\tuple{0,1,3}$, we use the frame rule with the post-call-site annotation $\tuple{0,0,2}$ to derive $\tuple{0,1,4(d-x)+11}$ as the pre-annotation:

$\{ \; \tuple{0, \;1,\; 4(d-x)+ 11} \; \}$  # $=\tuple{0, \; 1, \; 4(d-x)+9} \oplus \tuple{0, \; 0, \; 2}$

$\{ \;  \tuple{0, \; 1, \; 3} \;  \}$    # $=\tuple{0, \; 1, \; 1} \oplus \tuple{0, \; 0, \; 2}$

* ${\color{ACMDarkBlue} \{  \tuple{ 1,  2(d\!-\!x)\!+\!4 , 4(d\!-\!x)^2 \!+\! 22(d\!-\!x)\!+\!28  }  \} } ~ {\m{rdwalk}} ~{ \color{ACMDarkBlue} \{  \tuple{1,  0, 0}  \} }$: To analyze the recursive call with post-annotation $\tuple{1,1,1}$, we use the frame rule with the post-call-site annotation $\tuple{0,1,1}$ to derive $\tuple{1,2(d-x)+5,4(d-x)^2+26(d-x)+37}$ as the pre-annotation:

$\{ \; \tuple{1, 2(d\!-\!x)\!+\!5, 4(d\!-\!x)^2\!+\! 26(d\!-\!x)\!+\!37} \; \}$

# $=\tuple{1, 2(d\!-\!x)\!+\!4, 4(d\!-\!x)^2\!+\!22(d\!-\!x)\!+\!28} \oplus \tuple{0, 1,  4(d\!-\!x)\!+\!9}$

$\{ \;  \tuple{1, \; 1, \; 1} \;  \}$    # $=\tuple{1, \; 0, \; 0} \oplus \tuple{0, \; 1, \; 1}$

In <ref>, we present an automatic inference system for the expected-potential method that is extended with interval-valued bounds on higher moments, with support for moment-polymorphic recursion.

Soundness of the analysis

Unlike the classic potential method, the expected-potential method is not always sound when reasoning about the moments for cost accumulators in probabilistic programs.

func $\m{geo}$() begin  $\{ \; \tuple{1, \; 2^x } \;\}$
$\iassign{x}{\eadd{x}{1}}$;  $\{ \; \tuple{1, \; 2^{x-1} } \; \}$

# expected-potential method for lower bounds:

# $2^{x-1} < \sfrac{1}{2} \cdot (2^x+1) + \sfrac{1}{2} \cdot 0$

if $\kw{prob}(\sfrac{1}{2})$ then  $\{ \; \tuple{1, \; 2^{x} +1  } \;\}$
$\itick{1}$;  $\{ \; \tuple{1, \; 2^x} \; \}$

$\iinvoke{\m{geo}}$  $\{ \; \tuple{1, \; 0} \}$

A purely probabilistic loop with annotations for a lower bound on the first moment of the accumulated cost.

Consider the program in <ref> that describes a purely probabilistic loop that exits the loop with probability $\sfrac{1}{2}$ in each iteration.

The expected accumulated cost of the program should be one <cit.>.

However, the annotations in <ref> justify a potential function $2^x$ as a lower bound on the expected accumulated cost, no matter what value $x$ has at the beginning, which is apparently unsound.

Why does the expected-potential method fail in this case?

The short answer is that dualization only works for some problems:
upper-bounding the sum of nonnegative ticks is equivalent to
lower-bounding the sum of nonpositive ticks;
lower-bounding the sum of nonnegative ticks—the issue in
<ref>—is equivalent to upper-bounding the
sum of nonpositive ticks;
however, the two kinds of problems are inherently different <cit.>.

<cit.> presented a thorough study of the difficulty of lower-bound inference for probabilistic programs (which, also implies the difficulty of handling non-monotone costs, because upper-bounding the sum of nonnegative ticks is equivalent to lower-bounding the sum of nonpositive ticks).

Intuitively, the classic potential method for bounding the costs of non-probabilistic programs is a partial-correctness reasoning method, i.e., derived upper/lower bounds are sound if the analyzed program terminates <cit.>.

With probabilistic programs, many programs do not terminate definitely, but only almost surely, i.e., they terminate with probability one, but have some execution traces that are non-terminating.

The programs in <ref> are both almost-surely terminating.

For the expected-potential method, the potential at a program state can be seen as an average of potentials needed for all possible computations that continue from the state.

If the program state can lead to a non-terminating execution trace, the potential associated with that trace might be problematic, and as a consequence, the expected-potential method might fail.

Recent research <cit.> has employed the Optional Stopping Theorem (OST) from probability theory to address this soundness issue.

The classic OST provides a collection of sufficient conditions for reasoning about expected gain upon termination of stochastic processes, where the expected gain at any time is invariant.

By constructing a stochastic process for executions of probabilistic programs and setting the expected-potential function as the invariant, one can apply the OST to justify the soundness of the expected-potential function.

In a companion paper <cit.>, we study and propose an extension to the classic OST with a new sufficient condition that is suitable for reasoning about higher moments;

in this work, we prove the soundness of our central-moment inference for programs that satisfy this condition, and develop an algorithm to check this condition automatically (see <ref>).

\begin{align*}
S & \Coloneqq \iskip \mid \itick{c} \mid \iassign{x}{E} \mid \isample{x}{D} \mid \iinvoke{f} \mid \iloop{L}{S}  \\
& \mid \iprob{p}{S_1}{S_2} \mid \icond{L}{S_1}{S_2}  \mid S_1 ; S_2 \\
L & \Coloneqq \ctop \mid \cneg{L} \mid \cconj{L_1}{L_2} \mid \cle{E_1}{E_2} \\
E & \Coloneqq x \mid c \mid \eadd{E_1}{E_2} \mid \emul{E_1}{E_2} \\
D & \Coloneqq \kw{uniform}(a,b) \mid \cdots
\end{align*}

Syntax of , where $p \in [0,1]$, $a,b,c \in \bbR$, $a < b$, $x\in \mathsf{VID}$ is a variable, and $f \in \mathsf{FID}$ is a function identifier.

§ DERIVATION SYSTEM FOR HIGHER MOMENTS

In this section, we describe the inference system used by our analysis.

We first present a probabilistic programming language (<ref>).

We then introduce moment semirings to compose higher moments for a cost accumulator from two computations (<ref>).

We use moment semirings to develop our derivation system, which is presented as a declarative program logic (<ref>).

Finally, we sketch how we reduce the inference of a derivation to LP solving (<ref>).

§.§ A Probabilistic Programming Language

This paper uses an imperative arithmetic probabilistic programming language  that supports general recursion and continuous distributions, where program variables are real-valued.

We use the following notational conventions.

Natural numbers $\bbN$ exclude $0$, i.e., $\bbN \defeq \{ 1,2,3,\cdots\} \subseteq \bbZ^+ \defeq \{0,1,2,\cdots\}$.

The Iverson brackets $[\cdot]$ are defined by $[\varphi]=1$ if $\varphi$ is true and otherwise $[\varphi]=0$.

We denote updating an existing binding of $x$ in a finite map $f$ to $v$ by $f[x \mapsto v]$.

We will also use the following standard notions from probability theory: $\sigma$-algebras, measurable spaces, measurable functions, random variables, probability measures, and expectations.

<Ref> provides a review of those notions.
We include a review of those notions in the technical report <cit.>.

<ref> presents the syntax of , where the metavariables $S$, $L$, $E$, and $D$ stand for statements, conditions, expressions, and distributions, respectively.

Each distribution $D$ is associated with a probability measure $\mu_D \in \bbD(\bbR)$.

We write $\bbD(X)$ for the collection of all probability measures on the measurable space $X$.

For example, $\kw{uniform}(a,b)$ describes a uniform distribution on the interval $[a,b]$, and its corresponding probability measure is the integration of its density function $\mu_{\kw{uniform}(a,b)}(O) \defeq \int_O \frac{[a \le x \le b]}{b - a} dx $.

The statement “$\isample{x}{D}$” is a random-sampling assignment, which draws from the distribution $\mu_D$ to obtain a sample value and then assigns it to $x$.

The statement “$\iprob{p}{S_1}{S_2}$” is a probabilistic-branching statement, which executes $S_1$ with probability $p$, or $S_2$ with probability $(1-p)$.

The statement “$\iinvoke{f}$” makes a (possibly recursive) call to the function with identifier $f \in \mathsf{FID}$.

In this paper, we assume that the functions only manipulate states that consist of global program variables.

The statement $\itick{c}$, where $c \in \bbR$ is a constant, is used to define the cost model.

It adds $c$ to an anonymous global cost accumulator.

Note that our implementation supports local variables, function parameters, return statements, as well as accumulation of non-constant costs;

the restrictions imposed here are not essential, and are introduced solely to simplify the presentation.

We use a pair $\tuple{\scrD, S_{\mathsf{main}}}$ to represent an  program, where $\scrD$ is a finite map from function identifiers to their bodies and $S_{\mathsf{main}}$ is the body of the main function.

We present an operational semantics for  in <ref>.

§.§ Moment Semirings

As discussed in <ref>, we want to design a composition operation $\otimes$ and a combination operation $\oplus$ to compose and combine higher moments of accumulated costs such that

\begin{align*}
\phi(\sigma) & \!\sqsupseteq\! \expe_{\sigma' \!\sim\! \interp{S}(\sigma)} \![ \many{ \tuple{ C(\sigma,\sigma')^k}_{0\!\le\! k \!\le\! m}} \!\otimes\! \phi(\sigma')    ], \\
\phi_1(\sigma) \!\oplus\! \phi_2(\sigma) & \!\sqsupseteq\! \expe_{\sigma' \!\sim\! \interp{S}(\sigma)} \! [ \many{ \tuple{C(\sigma,\sigma')^k}_{0 \!\le\! k \!\le\! m} } \!\otimes\! (\phi_1(\sigma') \!\oplus\! \phi_2(\sigma')) ],
\end{align*}

where the expected-potential functions $\phi,\phi_1,\phi_2$ map program states to interval-valued vectors, $C(\sigma,\sigma')$ is the cost for the computation from $\sigma$ to $\sigma'$, and $m$ is the degree of the target moment.

In <ref>, we gave a definition of $\otimes$ and $\oplus$ suitable for first and second moments, respectively.

In this section, we generalize them to reason about upper and lower bounds of higher moments.

Our approach is inspired by the work of <cit.>, which develops a method to “lift” techniques for first moments to those for second moments.

Instead of restricting the elements of semirings to be vectors of numbers, we propose algebraic moment semirings that can also be instantiated with vectors of intervals, which we need for the interval-bound analysis that was demonstrated in <ref>.

The $m$-th order moment semiring $\calM_\calR^{(m)} = (|\calR|^{m+1},\oplus, \otimes, \azero, \aone)$ is parametrized by a partially ordered semiring $\calR = (|\calR|,{\le}, {+},{\cdot},0 ,1 )$,

\begin{align}
\many{\tuple{u_k}_{0 \le k \le m} } \!\oplus\! \many{\tuple{v_k}_{0 \le k \le m}} & \defeq \many{\tuple{u_k + v_k}_{0 \le k \le m}}, \\
\many{\tuple{u_k}_{0 \le k \le m}} \!\otimes\! \many{\tuple{v_k}_{0 \le k \le m}} & \defeq \many{ \tuple{ \textstyle \sum_{i=0}^k  \binom{k}{i} \!\times\! (u_i \cdot v_{k-i}) }_{0 \le k \le m} }, \label{Eq:MomentGenerating}
\end{align}

$\binom{k}{i}$ is the binomial coefficient; the scalar product $n \times u$ is an abbreviation for $\sum_{i=1}^n u$, for $n \in \bbZ^+, u \in \calR$; $\azero \defeq \tuple{0,0,\cdots,0}$; and $\aone \defeq \tuple{1 , 0 ,\cdots, 0 }$.

We define the partial order $\aord$ as the pointwise extension of the partial order $\le$ on $\calR$.

Intuitively, the definition of $\otimes$ in <ref> can be seen as the multiplication of two moment-generating functions for distributions with moments $\many{\tuple{u_k}_{0 \le k \le m}}$ and $\many{\tuple{v_k}_{0 \le k \le m}}$, respectively.

We prove a composition property for moment semirings.

For all $u,v \in \calR$, it holds that

\[
\many{\tuple{(u+v)^k}_{0 \le k \le m}} = \many{\tuple{u^k}_{0 \le k \le m}} \otimes \many{\tuple{v^k}_{0 \le k \le m}},
\]

where $u^n$ is an abbreviation for $\prod_{i=1}^n u$, for $n \in \bbZ^+, u \in \calR$.

§.§ Inference Rules

We present the derivation system as a declarative program logic that uses moment semirings to enable compositional reasoning and moment-polymorphic recursion.

Interval-valued moment semirings

Our derivation system infers upper and lower bounds simultaneously, rather than separately,
which is essential for non-monotone costs.

Consider a program “$\itick{-1} ; S$” and suppose that we have $\tuple{1,2,5}$ and $\tuple{1,{-2},5}$ as the upper and lower bound on the first two moments of the cost for $S$, respectively.

If we only use the upper bound, we derive $\tuple{1,{-1},1} \otimes \tuple{1,2,5} = \tuple{1,1,2}$, which is not an upper bound on the moments of the cost for the program;
if the actual moments of the cost for $S$ are $\tuple{1,0,5}$,
then the actual moments of the cost for “$\itick{-1}; S$” are $\tuple{1,{-1},1} \otimes \tuple{1,0,5} = \tuple{1,{-1}, 4} \not\le \tuple{1,1,2}$.

Thus, in the analysis, we instantiate moment semirings with the interval domain $\calI$.

For the program “$\itick{-1}; S$”, its interval-valued bound on the first two moments is
\tuple{[1,1], [{-1},{-1}],[1,1]} \otimes \tuple{[1,1], [{-2},2], [5,5]}
=  \tuple{[1,1], [{-3},1], [ 2, 10]}.

Template-based expected-potential functions

The basic approach to automated inference using potential functions is to introduce a template for the expected-potential functions.

Let us fix $m \in \bbN$ as the degree of the target moment.

Because we use $\calM_\calI^{(m)}$-valued expected-potential functions whose range is vectors of intervals, the templates are vectors of intervals whose ends are represented symbolically.

In this paper, we represent the ends of intervals by polynomials in $\bbR[\mathsf{VID}]$ over program variables.

More formally, we lift the interval semiring $\calI$ to a symbolic interval semiring $\calP\calI$ by representing the ends of the $k$-th interval by polynomials in $\bbR_{kd}[\mathsf{VID}]$ up to degree $kd$ for some fixed $d \in \bbN$.

Let $\calM_{\calP\calI}^{(m)}$ be the $m$-th order moment semiring instantiated with the symbolic interval semiring.

Then the potential annotation is represented as $Q = \many{\tuple{[L_k,U_k]}_{0 \le k \le m}} \in \calM_{\calP\calI}^{(m)}$, where $L_k$'s and $U_k$'s are polynomials in $\bbR_{kd}[\mathsf{VID}]$.

$Q$ defines an $\calM_\calI^{(m)}$-valued expected-potential function $\phi_Q(\sigma) \defeq \many{ \tuple{[ \sigma(L_k),\sigma(U_k) ]}_{0 \le k \le m} }$, where $\sigma$ is a program state, and $\sigma(L_k)$ and $\sigma(U_k)$ are $L_k$ and $U_k$ evaluated over $\sigma$, respectively.

Inference rules

We formalize our derivation system for moment analysis in a Hoare-logic style.

The judgment has the form $\Delta \vdash_h \{\Gamma;Q\}~S~\{\Gamma';Q'\}$, where
$S$ is a statement, $\{\Gamma;Q\}$ is a precondition, $\{\Gamma';Q'\}$ is a postcondition,
$\Delta = \tuple{\Delta_k}_{0 \le k \le m}$ is a context of function specifications,
and $h \in \bbZ^+$ specifies some restrictions put on $Q,Q'$ that we will explain later.

The logical context $\Gamma\!:\!(\mathsf{VID} \!\to\! \bbR) \!\to\! \set{\top,\bot}$ is a predicateon program states that describes reachable states at a program point.

The potential annotation $Q \in \calM_{\calP\calI}^{(m)}$ specifies a map from program states to the moment semiring that is used to define interval-valued expected-potential functions.

The semantics of the triple $\{\cdot;Q\}~S~\{\cdot;Q'\}$ is that if the rest of the computation after executing $S$ has its moments of the accumulated cost bounded by $\phi_{Q'}$, then the whole computation has its moments of the accumulated cost bounded by $\phi_Q$.

The parameter $h$ restricts all $i$-th-moment components in $Q,Q'$, such that $i < h$, to be $[0,0]$.

We call such potential annotations $h$-restricted;

this construction is motivated by an observation from <ref>, where we illustrated the benefits of carrying out interprocedural analysis using an “elimination sequence” of annotations for recursive function calls, where the successive annotations have a greater number of zeros, filling from the left.

Function specifications are valid pairs of pre- and post-conditions for all declared functions in a program.

For each $k$, such that $0 \!\le\! k \!\le\! m$, and each function $f$,
a valid specification $(\Gamma;Q,\Gamma';Q') \in \Delta_k(f)$ is justified by the judgment $\Delta \vdash_k \{\Gamma;Q\}~\scrD(f)~\{\Gamma';Q'\}$, where $\scrD(f)$ is the function body of $f$, and $Q,Q'$ are $k$-restricted.

To perform context-sensitiveinterprocedural analysis, a function can have multiple specifications.

Q-Tick    Q = [c^k,c^k]_0 ≤k ≤m ⊗Q'  Δ⊢_h { Γ;Q } c { Γ;Q' }

Q-Sample Γ=x ∈supp(μ_D) Γ'

Q = _x ∼μ_D[Q']  Δ⊢_h {  Γ;Q   } xD { Γ' ;Q' }

Q-Loop Δ⊢_h { Γ∧L;Q } S { Γ;Q }   Δ⊢_h { Γ;Q }  LS  { Γ∧ L; Q }
(Γ;Q , Γ';Q') ∈Δ_m(f)

Δ⊢_m { Γ; Q }  f  { Γ'; Q' }

Q-Seq Δ⊢_h {Γ;Q } S_1 {Γ';Q' }

Δ⊢_h { Γ';Q' } S_2 {Γ”;Q”}   Δ⊢_h { Γ;Q } S_1;S_2 {Γ”;Q”}
Q-Call-Poly h < m
Δ_h(f) = (Γ;Q_1,Γ';Q_1')

Δ⊢_h+1 { Γ;Q_2} (f) { Γ';Q_2' }    Δ⊢_h { Γ; Q_1 ⊕Q_2 } f { Γ'; Q_1' ⊕Q_2' }

Q-Prob  Δ⊢_h { Γ;Q_1 } S_1 { Γ';Q' }
Δ⊢_h { Γ;Q_2 }  S_2 { Γ';Q' }
Q = P ⊕R
P = [p,p],[0,0],⋯,[0,0] ⊗Q_1
R = [1-p,1-p],[0,0],⋯,[0,0] ⊗Q_2  Δ⊢_h { Γ;Q }  pS_1S_2 { Γ'; Q' }

Q-Weaken Δ⊢_h { Γ_0;Q_0 } S { Γ_0';Q_0' }  ΓΓ_0  Γ_0' Γ'  ΓQ ⊒Q_0  Γ_0' Q_0' ⊒Q'   Δ⊢_h {Γ;Q } S {Γ';Q'}

Selected inference rules of the derivation system.

<ref> presents some of the inference rules.

The rule (Q-Tick) is the only rule that deals with costs in a program.

To accumulate the moments of the cost, we use the $\otimes$ operation in the moment semiring $\calM_{\calP\calI}^{(m)}$.

The rule (Q-Sample) accounts for sampling statements.

Because “$x \!\sim\! D$” randomly assigns a value to $x$ in the support of distribution $D$, we quantify $x$ out universally from the logical context.

To compute $Q = \expe_{x \sim \mu_D}[Q']$, where $x$ is drawn from distribution $D$, we assume the moments for $D$ are well-defined and computable,
and substitute $x^i$, $i \in \bbN$ with the corresponding moments in $Q'$.

We make this assumption because every component of $Q'$ is a polynomial over program variables.

For example, if $D = \kw{uniform}({-1},2)$, we know the following facts

\[
\expe_{x \sim \mu_D}[x^0] \!=\! 1, \expe_{x \sim \mu_D}[x^1] \!=\! \sfrac{1}{2},\expe_{x \sim \mu_D}[x^2] \!=\! 1,\expe_{x \sim \mu_D}[x^3] \!=\! \sfrac{5}{4}.
\]

Then for $Q' = \langle [1,1], [1+x^2,xy^2+x^3y] \rangle$, by the linearity of expectations, we compute $Q=\expe_{x \sim \mu_D}[Q']$ as follows:

\begin{aligned}
\expe_{x \sim \mu_D}[Q'] & = \langle [1,1], [\expe_{x \sim \mu_D}[1+x^2], \expe_{x \sim \mu_D}[xy^2+x^3y] \rangle \\[-3pt]
& = \langle [1,1], [1 + \expe_{x \sim \mu_D}[x^2], y^2 \expe_{x \sim \mu_D}[x] + y \expe_{x \sim \mu_D}[x^3]] \rangle \\[-3pt]
& = \langle [1,1], [2, \sfrac{1}{2} \cdot  y^2+\sfrac{5}{4} \cdot y] \rangle.
\end{aligned}

The other probabilistic rule (Q-Prob) deals with probabilistic branching.

Intuitively, if the moments of the
execution of $S_1$ and $S_2$ are $q_1$ and $q_2$, respectively, and those of the accumulated cost of the computation after the branch statement is bounded by $\phi_{Q'}$,
then the moments for the whole computation should be bounded by a “weighted average” of $(q_1 \otimes \phi_{Q'})$ and $(q_2 \otimes \phi_{Q'})$, with respect to the branching probability $p$.

We implement the weighted average by the combination operator $\oplus$ applied to
\tuple{ [p,p],[0,0],\cdots,[0,0] } \otimes q_1 \otimes \phi_{Q'}
$\tuple{ [1-p,1-p],[0,0],\cdots,[0,0] } \otimes q_2 \otimes \phi_{Q'}$,
because the $0$-th moments denote probabilities.

The rules (Q-Call-Poly) and (Q-Call-Mono) handle function calls.

Recall that in <ref>, we use the $\oplus$ operator to combine multiple potential functions for a function to reason about recursive function calls.

The restriction parameter $h$ is used to ensure that the derivation system only needs to reason about finitely many post-annotations for each call site.

In rule (Q-Call-Poly), where $h$ is smaller than the target moment $m$,
we fetch the pre- and post-condition $Q_1,Q_1'$ for the function $f$ from the specification context $\Delta_h$.

We then combine it with a frame of $(h+1)$-restricted potential annotations $Q_2,Q_2'$ for the function $f$.

The frame is used to account for the interval bounds on the moments for the computation after the function call for most non-tail-recursive programs.

When $h$ reaches the target moment $m$, we use the rule (Q-Call-Mono) to reason moment-monomorphically, because setting $h$ to $m+1$ implies that the frame can only be $\tuple{[0,0],[0,0],\cdots,[0,0]}$.

The structural rule (Q-Weaken) is used to strengthen the pre-condition and relax the post-condition.

The entailment relation $\Gamma \models \Gamma'$ states that the logical implication $\Gamma \implies \Gamma'$ is valid.

In terms of the bounds on higher moments for cost accumulators, if the triple $\{\cdot;Q\}~S~\{\cdot;Q'\}$ is valid, then we can safely widen the intervals in the pre-condition $Q$ and narrow the intervals in the post-condition $Q'$.

func $\mathsf{rdwalk}$() begin
${\color{ACMDarkBlue} \{ \; x \!<\! d \!+\! 2; \; \langle [1,1], \; [2(d \!-\! x), 2(d\!-\!x)\!+\!4], }$

${\color{ACMDarkBlue} \enskip [4(d \!-\! x)^2 + 6(d \!-\! x) \!-\! 4, 4(d\!-\!x)^2\!+\!22(d\!-\!x)\!+\!28] \rangle  \; \} }$

if $\cbin{<}{x}{d}$ then
${\color{ACMDarkBlue} \{ \;  x  \!<\!  d;  \; \langle [1,1], \; [2(d \!-\! x), 2(d\!-\!x)\!+\!4], }$

${\color{ACMDarkBlue} \enskip [4(d \!-\! x)^2 \!+\! 6(d\!-\! x) \!-\! 4, 4(d\!-\!x)^2\!+\!22(d\!-\!x)\!+\!28 ] \rangle \; \}  }$

$\isample{t}{ \kw{uniform}({-1},2)}$;

${\color{ACMDarkBlue} \{ \; x  \!<\!  d  \!\wedge\!  t  \!\le\!  2 ; \; \langle [1,1], \; [2(d \!-\! x \!-\! t) \!+\! 1, 2(d\!-\!x\!-\!t)\!+\!5], }$

${\color{ACMDarkBlue} \enskip [4(d \!-\! x \!-\! t)^2 + 10(d \!-\! x \!-\! t) \!-\! 3, 4(d\!-\!x\!-\!t)^2\!+\!26(d\!-\!x\!-\!t)\!+\!37) ] \rangle \;  \} }$

$ \iassign{x}{  \eadd{x}{t} }$;

${\color{ACMDarkBlue} \{ \; x  \!<\!  d  \!+\!  2; \; \langle [1,1], \; [2(d \!-\! x) \!+\! 1, 2(d\!-\!x)\!+\!5], }$

${\color{ACMDarkBlue} \enskip [4(d \!-\! x)^2 \!+\! 10(d \!-\! x) \!-\! 3, 4(d\!-\!x)^2\!+\!26(d\!-\!x)\!+\!37 ] \rangle \; \} }$

call $\mathsf{rdwalk}$;

${\color{ACMDarkBlue} \{ \; \top; \; \tuple{ [1,1], \; [1,1], \; [1,1] } \; \} }$

${\color{ACMDarkBlue} \{ \; \top; \; \tuple{[1,1],\;[0,0],\;[0,0]} \; \} } $

The $\mathsf{rdwalk}$ function with annotations for the interval-bounds on the first and second moments.

<ref> presents the logical context and the complete potential annotation for the first and second moments for the cost accumulator $\id{tick}$ of the $\mathsf{rdwalk}$ function from <ref>.

Similar to the reasoning in <ref>, we can justify the derivation using moment-polymorphic recursion and the moment bounds for $\m{rdwalk}$ with post-annotations
$\tuple{[0,0],[1,1],[1,1]}$ and $\tuple{[0,0],[0,0],[2,2]}$.

§.§ Automatic Linear-Constraint Generation

We adapt existing techniques <cit.> to automate our inference system by
(i) using an abstract interpreter to infer logical contexts,
(ii) generating templates and linear constraints by inductively applying the derivation rules to the analyzed program, and
(iii) employing an off-the-shelf LP solver to discharge the linear constraints.

During the generation phase, the coefficients of monomials in the polynomials from the ends of the intervals in every qualitative context $Q \in \calM_{\calP\calI}^{(m)}$ are recorded as symbolic names, and the inequalities among those coefficients—derived from the inference rules in <ref>—are emitted to the LP solver.

Generating linear constraints

<ref> demonstrates the generation process for some of the bounds in <ref>.

Let $B_k$ be a vector of monomials over program variables $\mathsf{VID}$ of degree up to $k$.
Then a polynomial $\sum_{b \in B_k} q_b \cdot b$, where $q_b \in \bbR$ for all $b \in B_k$, can be represented as a vector of its coefficients $(q_b)_{b \in B_k}$.

We denote coefficient vectors by uppercase letters, while we use lowercase letters as names of the coefficients.

We also assume that the degree of the polynomials for the $k$-th moments is up to $k$.

p^tk_1 = u^tk_1

q^tk_1 = u^tk_1 + v^tk_1

q^tk_x = v^tk_x

q^tk_N = v^tk_N

q^tk_r = v^tk_r

t^tk_1 = w^tk_1 + 2v^tk_1 + u^tk_1

t^tk_x = w^tk_x + 2v^tk_x

t^tk_x^2 = w^tk_x^2

Δ⊢{ ⊤; (P^tk,Q^tk,T^tk) } 1 {⊤; (U^tk,V^tk,W^tk) }
p^sa_1 = u^sa_1

q^sa_1 = v^sa_1 + 12 ·v^sa_r

q^sa_x = v^sa_x

q^sa_N = v^sa_N

q^sa_r = 0

t^sa_1 = w^sa_1 + 12 ·w^sa_r + 1 ·w^sa_r^2

t^sa_x = w^sa_x + 12 ·w^sa_r ·x

t^sa_x^2 = w^sa_x^2

Δ⊢{ x<N; (P^sa,Q^sa,T^sa) } runiform(-1,2)  {x < N ∧r ≤2; (U^sa,V^sa,W^sa) }

Generate linear constraints, guided by inference rules.

For (Q-Tick), we generate constraints that correspond to the composition operation $\otimes$ of the moment semiring.

For example, the second-moment component should satisfy

\begin{align*}
\sum_{b \in B_2} t^{tk}_b \cdot b & = \text{the second-moment component of}~ \big( (1,1,1) \\
& \otimes (\sum_{b \in B_0} u^{tk}_b \cdot b,\sum_{b \in B_1} v^{tk}_b \cdot b,\sum_{b \in B_2} w^{tk}_b \cdot b) \big) \\
& = \sum_{b \in B_2} w^{tk}_b \cdot b + 2 \cdot \sum_{b \in B_1} v^{tk}_b \cdot b + \sum_{b \in B_0} u^{tk}_b \cdot b.
\end{align*}

Then we extract $t^{tk}_1 = w^{tk}_1 + 2v^{tk}_1 + u^{tk}_1$ for $b=1$, and $t^{tk}_x = w^{tk}_x + 2v^{tk}_x$ for $b=x$, etc.

For (Q-Sample), we generate constraints to perform “partial evaluation” on the polynomials by substituting $r$ with the moments of $\kw{uniform}({-1},2)$.

As we discussed in <ref>, let $D$ denote $\kw{uniform}({-1},2)$, then $\expe_{r\sim D}[w^{sa}_r \cdot r] = w^{sa}_r \cdot \expe_{r \sim D}[r] = \sfrac{1}{2} \cdot w^{sa}_r$, $\expe_{r \sim D}[w^{sa}_{r^2} \cdot r^2] = w^{sa}_{r^2} \cdot \expe_{r \sim D}[r^2] = 1 \cdot w^{sa}_{r^2}$.

Then we generate a constraint $t^{sa}_1 = w^{sa}_1 + \sfrac{1}{2} \cdot w^{sa}_r + 1 \cdot w^{sa}_{r^2}$ for $t^{sa}_1$.

The loop rule (Q-Loop) involves constructing loop invariants $Q$, which is in general a non-trivial problem for automated static analysis.

Instead of computing the loop invariant $Q$ explicitly, our system represents $Q$ directly as a template with unknown coefficients,
then uses $Q$ as the post-annotation to analyze the loop body and obtain a pre-annotation,
and finally generates linear constraints that indicate the pre-annotation equals to $Q$.

The structural rule (Q-Weaken) can be applied at any point during the derivation.
In our implementation, we apply it where the control flow has a branch, because different branches might have different costs.

To handle the judgment $\Gamma \models Q \sqsupseteq Q'$, i.e., to generate constraints that ensure one interval is always contained in another interval, where the ends of the intervals are polynomials,
we adapt the idea of rewrite functions <cit.>.

Intuitively, to ensure that $[L_1,U_1] \sqsupseteq [L_2,U_2]$, i.e., $L_1 \le L_2$ and $U_2 \le U_1$, under the logical context $\Gamma$, we generate constraints indicating that there exist two polynomials $T_1,T_2$ that are always nonnegative under $\Gamma$, such that $L_1=L_2+T_1$ and $U_1=U_2-T_2$.

Here, $T_1$ and $T_2$ are like slack variables, except that because all quantities are polynomials, they are too (i.e., slack polynomials).

In our implementation, $\Gamma$ is a set of linear constraints over program variables of the form $\calE \ge 0$ , then we can represent $T_1,T_2$ by conical combinations (i.e., linear combinations with nonnegative scalars) of expressions $\calE$ in $\Gamma$.

We demonstrate linear-constraint generation for the upper bound on the first moment for the sampling statement $\isample{x}{\kw{uniform}({-1},2)}$ with a pre-annotation $\tuple{[0,0], [0,q_{x^2} \cdot x^2 + q_{x} \cdot x + q_{1} \cdot 1]}$ and a post-annotation $\tuple{[0,0], [0,q_{x^2}' \cdot x^2 + q_{x}' \cdot x + q_{1}' \cdot 1]}$,

where we use polynomials of $x$ up to degree 2 as the templates, and $q_{x^2},q_x,q_1,q_{x^2}',q_x',q_1'$ are unknown numeric coefficients.

By (Q-Sample), we generate constraints to perform “partial evaluation” on the polynomials by substituting $x$ with the moments of $\kw{uniform}({-1},2)$.

Let $D$ denote $\kw{uniform}({-1},2)$. Because

\[
\expe_{x \sim \mu_D}[q_{x^2}' \cdot x^2 + q_{x}' \cdot x + q_{1}' \cdot 1] = (q_{x^2}' \cdot 1 + q_{x}' \cdot \sfrac{1}{2} + q_{1}' \cdot 1),
\]

we generate these linear constraints:

\[
q_{x^2}=0, \quad q_x=0, \quad q_1 = q'_{x^2} + q'_x \cdot \sfrac{1}{2} + q_1'.
\]

Constraint generation for other inference rules is similar to the process we describe in <ref>.

For example, let us consider the loop rule (Q-Loop).

Instead of computing the loop invariant $Q$ explicitly, our system represents $Q$ directly as a template with unknown coefficients,
then uses $Q$ as the post-annotation to analyze the loop body and obtain a pre-annotation,
and finally generates linear constraints that indicate that the pre-annotation equals $Q$.

Details of the linear-constraint generation of our system are included in the technical report <cit.>.

Solving linear constraints

The LP solver not only finds assignments to the coefficients that satisfy the constraints, it can also optimize a linear objective function.

In our central-moment analysis, we construct an objective function that tries to minimize imprecision.

For example, let us consider upper bounds on the variance.
We randomly pick a concrete valuation of program variables that satisfies the pre-condition (e.g., $d>0$ in <ref>), and then substitute program variables with the concrete valuation in the polynomial for the upper bound on the variance (obtained from bounds on the raw moments).

The resulting linear combination of coefficients, which we set as the objective function, stands for the variance under the concrete valuation.

Thus, minimizing the objective function produces the most precise upper bound on the variance under the specific concrete valuation.

Also, we can extract a symbolic upper bound on the variance using the assignments to the coefficients.

Because the derivation of the bounds only uses the given pre-condition, the symbolic bounds apply to all valuationsof program variables that satisfy the pre-condition.

§ SOUNDNESS OF HIGHER-MOMENT ANALYSIS

In this section, we study the soundness of our derivation system for higher-moment analysis.

We first present a Markov-chain semantics for the probabilistic programming language  to reason about how stepwise costs contribute to the global accumulated cost (<ref>).

We then formulate higher-moment analysis with respect to the semantics and prove the soundness of our derivation system for higher-moment analysis based on a recent extension to the Optional Stopping Theorem (<ref>).

Finally, we sketch the algorithmic approach for ensuring the soundness of our analysis (<ref>).

§.§ A Markov-Chain Semantics

Operational semantics

We start with a small-step operational semantics with continuations, which we will use later to construct the Markov-chain semantics.

We follow a distribution-based approach <cit.>
to define an operational cost semantics for .

Full details of the semantics are included in
the technical report <cit.>.

A program configuration $\sigma \in \Sigma$ is a quadruple $\tuple{\gamma,S,K,\alpha}$ where $\gamma : \mathsf{VID} \to \bbR$ is a program state that maps variables to values, $S$ is the statement being executed, $K$ is a continuation that describes what remains to be done after the execution of $S$, and $\alpha \in \bbR$ is the global cost accumulator.

An execution of an  program $\tuple{\scrD,S_\mathsf{main}}$ is initialized with $\tuple{\lambda\_. 0, S_{\mathsf{main}}, \kstop,0}$, and the termination configurations have the form $\tuple{\_,\iskip,\kstop,\_}$,
where $\kstop$ is an empty continuation.

Different from a standard semantics where each program configuration steps to at most one new configuration, a probabilistic semantics may pick several different new configurations.

The evaluation relation for  has the form $\sigma \mapsto \mu$ where $\mu \in \bbD(\Sigma)$ is a probability measure over configurations.

Below are two example rules.

The rule (E-Prob) constructs a distribution whose support has exactly two elements, which stand for the two branches of the probabilistic choice.

We write $\delta(\sigma)$ for the Dirac measure at $\sigma$, defined as $\lambda A.[\sigma\in A]$ where $A$ is a measurable subset of $\Sigma$.

We also write $p \cdot \mu_1 + (1-p)\cdot \mu_2$ for a convex combination of measures $\mu_1$ and $\mu_2$ where $p \in [0,1]$, defined as $\lambda A. p \cdot \mu_1(A) + (1-p) \cdot \mu_2(A)$.

The rule (E-Sample) “pushes” the probability distribution of $D$ to a distribution over post-sampling program configurations.

*[right=(E-Prob)]S=pS_1S_2  γ, S, K,α ↦p ·δ(γ,S_1,K,α) + (1-p) ·δ(γ,S_2,K,α)
*[right=(E-Sample)]  γ, x∼D, K,α ↦λA. μ_D({ r |γ[x ↦ r],,K,α ∈A})

Suppose that a random sampling statement is being executed, i.e., the current configuration is

$\tuple{\{ t \mapsto t_0 \}, (\isample{t}{\kw{uniform}(-1,2)}),K_0,\alpha_0}.$

The probability measure for the uniform distribution is $\lambda O. \int_O \frac{[-1 \le x \le 2]}{3}   dx$.

Thus, by the rule (E-Sample), we derive the post-sampling probability measure over configurations:

\lambda A. \int_{\bbR} [\tuple{ \{ t \mapsto r \}, \iskip, K_0,\alpha_0  } \in A] \cdot \frac{[-1 \le r \le 2]}{3} dr.

A Markov-chain semantics

In this work, we harness Markov-chain-based reasoning <cit.> to develop a Markov-chain cost semantics for , based on the evaluation relation $\sigma \mapsto \mu$.

An advantage of this approach is that it allows us to study how the cost of every single evaluation step contributes to the accumulated cost at the exit of the program.

Details of this semantics are included in
the technical report <cit.>.

Let $(\Omega,\calF,\prob)$ be the probability space where $\Omega \defeq \Sigma^{\bbZ^+}$ is the set of all infinite traces over program configurations, $\calF$ is a $\sigma$-algebra on $\Omega$, and $\prob$ is a probability measure on $(\Omega,\calF)$ obtained by the evaluation relation $\sigma \mapsto \mu$ and the initial configuration $\tuple{\lambda\_.0, S_{\mathsf{main}},\kstop,0}$.

Intuitively, $\prob$ specifies the probability distribution over all possible executions of a probabilistic program.

The probability of an assertion $\theta$ with respect to $\prob$, written $\prob[\theta]$, is defined as $\prob(\{ \omega \mid \theta(\omega)~\text{is true} \})$.

To formulate the accumulated cost at the exit of the program, we define the stopping time $T:\Omega \to \bbZ^+ \cup \{\infty\}$ of a probabilistic program as a random variable on the probability space $(\Omega,\calF,\prob)$ of program traces:

\[
T(\omega) \defeq \inf \{ n \in \bbZ^+ \mid \omega_n = \tuple{\_,\iskip,\kstop,\_} \},
\]

i.e., $T(\omega)$ is the number of evaluation steps before the trace $\omega$ reaches some termination configuration $\tuple{\_,\iskip,\kstop,\_}$.

We define the accumulated cost $A_T : \Omega \to \bbR$ with respect to the stopping time $T$ as

\[
A_T(\omega) \defeq A_{T(\omega)}(\omega),
\]

where $A_n : \Omega \to \bbR$ captures the accumulated cost at the $n$-th evaluation step for $n \in \bbZ^+$, which is defined as

\[
A_n(\omega) \defeq \alpha_n ~\text{where}~\omega_n =\tuple{\_,\_,\_,\alpha_n}.
\]

The $m$-th moment of the accumulated cost is given by the expectation $\expe[A_T^m]$ with respect to $\prob$.

§.§ Soundness of the Derivation System

Proofs for this section are included in
the technical report <cit.>.

The expected-potential method for moment analysis

We fix a degree $m \in \bbN$ and let $\calM_\calI^{(m)}$ be the $m$-th order moment semiring instantiated with the interval semiring $\calI$.

We now define $\calM_\calI^{(m)}$-valued expected-potential functions.

A measurable map $\phi : \Sigma \to \calM_\calI^{(m)}$ is said to be an expected-potential function if

* $\phi(\sigma) = \aone$ if $\sigma=\tuple{\_,\iskip,\kstop,\_}$, and

* $\phi(\sigma) \sqsupseteq \expe_{\sigma' \sim {\mapsto}(\sigma)}[\many{\big([(\alpha'-\alpha)^k,(\alpha'-\alpha)^k]\big)_{0 \le k \le m}} \otimes \phi(\sigma')]$ where $\sigma=\tuple{\_,\_,\_,\alpha},\sigma'=\tuple{\_,\_,\_,\alpha'}$ for all $\sigma \in \Sigma$.

Intuitively, $\phi(\sigma)$ is an interval bound on the moments for the accumulated cost of the computation that continues from the configuration $\sigma$.

We define $\mathbf{\Phi}_n$ and $\mathbf{Y}_n$, where $n \in \bbZ^+$, to be $\calM_\calI^{(m)}$-valued random variables on the probability space $(\Omega,\calF,\prob)$ of the Markov-chain semantics as

\[
\mathbf{\Phi}_n(\omega) \defeq \phi(\omega_n),
\mathbf{Y}_n(\omega) \defeq \many{ \tuple{ [A_n(\omega)^k, A_n(\omega)^k] }_{0 \le k \le m} } \otimes \mathbf{\Phi}_n(\omega)  .
\]

In the definition of $\mathbf{Y}_n$, we use $\otimes$ to compose the powers of the accumulated cost at step $n$ and the expected potential function that stands for the moments of the accumulated cost for the rest of the computation.

By the properties of potential functions, we can prove that $\expe[\mathbf{Y}_{n+1} \mid \mathbf{Y}_n] \sqsubseteq \mathbf{Y}_n$ almost surely, for all $n \in \bbZ^+$.

We call $\{\mathbf{Y}_n\}_{n \in \bbZ^+}$ a moment invariant.

Our goal is to establish that $\expe[\mathbf{Y}_T] \sqsubseteq \expe[\mathbf{Y}_0]$, i.e., the initial interval-valued potential $\expe[\mathbf{Y}_0] = \expe[\aone \otimes \mathbf{\Phi}_0] = \expe[\mathbf{\Phi}_0]$ brackets the higher moments of the accumulated cost $\expe[\mathbf{Y}_T] = \expe[\many{ \tuple{ [A_T^k, A_T^k] }_{0 \le k \le m}  } \otimes \aone] = \many{\tuple{ [\expe[A_T^k], \expe[A_T^k] ]  }_{0 \le k \le m} }$.

The soundness of the derivation system is proved with respect to the Markov-chain semantics.

Let $\norm{\many{\tuple{[a_k,b_k]}_{0 \le k \le m} }}_\infty \defeq \max_{0 \le k \le m} \{ \max\{ |a_k|,|b_k| \} \} $.

Let $\tuple{\scrD,S_\mathsf{main}}$ be a probabilistic program.
Suppose $\Delta \vdash \{\Gamma;Q\}~S_{\mathsf{main}}~\{\Gamma';\aone\}$, where $Q \in \calM_{\calP\calI}^{(m)}$ and the ends of the $k$-th interval in $Q$ are polynomials in $\bbR_{kd}[\mathsf{VID}]$.

Let $\{\mathbf{Y}_n\}_{n \in \bbZ^+}$ be the moment invariant extracted from the Markov-chain semantics with respect to the derivation of $\Delta \vdash \{\Gamma;Q\}~S_{\mathsf{main}}~\{\Gamma';\aone\}$.
If the following conditions hold:

*  $\expe[T^{md}] < \infty$, and

*  there exists $C \ge 0$ such that for all $n \in \bbZ^+$, $\norm{\mathbf{Y}_n}_\infty \le C \cdot (n+1)^{md}$ almost surely,

$\many{\tuple{[\expe[A_T^k], \expe[A_T^k] ]}_{0 \le k \le m}} \aord \phi_Q(\lambda\_.0)$.

The intuitive meaning of $\many{\tuple{[\expe[A_T^k], \expe[A_T^k]]}_{0 \le k \le m}} \aord \phi_Q(\lambda\_.0)$ is that the moment $\expe[A_T^k]$ of the accumulated cost upon program termination is bounded by the interval in the $k$th-moment component of $\phi_Q(\lambda\_.0)$, where $Q$ is the quantitative context and $\lambda\_.0$ is the initial state.

As we discussed in <ref> and <ref>, the expected-potential method is not always sound for deriving bounds on higher moments for cost accumulators in probabilistic programs.

The extra conditions <ref><ref> and <ref> impose
constraints on the analyzed program and the expected-potential function, which allow us to reduce the soundness to the optional stopping problem from probability theory.

Optional stopping

Let us represent the moment invariant $\{\mathbf{Y}_n\}_{n \in \bbZ^+}$ as

\[
\{ \tuple{[L^{(0)}_n, U^{(0)}_n], [L^{(1)}_n, U^{(1)}_n], \cdots, [L^{(m)}_n, U^{(m)}_n]  }  \}_{n \in \bbZ^+},
\]

where $L^{(k)}_n, U^{(k)}_n : \Omega \to \bbR$ are real-valued random variables on the probability space $(\Omega,\calF,\prob)$ of the Markov-chain semantics, for $n \in \bbZ^+$, $0 \le k \le m$.

We then have the observations below as direct corollaries of <ref>:

* For any $k$, the sequence $\{U_n^{(k)}\}_{n \in \bbZ^+}$ satisfies $\expe[U_{n+1}^{(k)} \mid U_n^{(k)}] \le U_n^{(k)}$ almost surely, for all $n \in \bbZ^+$, and we want to find sufficient conditions for $\expe[U_T^{(k)}] \le \expe[U_0^{(k)}]$.

* For any $k$, the sequence $\{L_n^{(k)}\}_{n \in \bbZ^+}$ satisfies $\expe[L_{n+1}^{(k)} \mid L_n^{(k)}] \ge L_n^{(k)}$ almost surely, for all $n \in \bbZ^+$, and we want to find sufficient conditions for $\expe[L_T^{(k)}] \ge \expe[L_0^{(k)}]$.

These kinds of questions can be reduced to optional stopping problem from probability theory.

Recent research <cit.> has used the Optional Stopping Theorem (OST) from probability theory to establish sufficient conditions for the soundness for analysis of probabilistic programs.

However, the classic OST turns out to be not suitable for higher-moment analysis.

We extend OST with a new sufficient condition that allows us to prove <ref>.
We discuss the details of our extended OST in a companion paper <cit.>;
in this work, we focus on the derivation system for central moments.

We also include a brief discussion on OST in <Ref>.

§.§ An Algorithm for Checking Soundness Criteria

Termination Analysis

We reuse our system for automatically deriving higher moments, which we developed in <ref>, for checking if $\expe[T^{md}] < \infty$ (<ref><ref>).

To reason about termination time, we assume that every program statement increments the cost accumulator by one.

For example, the inference rule
(Q-Sample) becomes

Γ= x ∈supp(μ_D) Γ'

Q = 1,1,⋯,1 ⊗_x ∼μ_D[Q']

Δ⊢{Γ;Q} xD {Γ';Q'}

However, we cannot apply <ref> for the soundness of the termination-time analysis, because that would introduce a circular dependence.

Instead, we use a different proof technique to reason about $\expe[T^{md}]$, taking into account the monotonicity of the runtimes.

Because upper-bound analysis of higher moments of runtimes has been studied by <cit.>, we skip the details, but include them in
the technical report <cit.>.

Boundedness of $\norm{Y_n}_\infty$, $n \in \bbZ^+$

To ensure that the condition in <ref><ref> holds, we check if the analyzed program satisfies the bounded-update property: every (deterministic or probabilistic) assignment to a program variable updates the variable with a change bounded by a constant $C$ almost surely.

Then the absolute value of every program variable at evaluation step $n$ can be bounded by $C \cdot n = O(n)$.
Thus, a polynomial up to degree $\ell \in \bbN$ over program variables can be bounded by $O(n^\ell)$ at evaluation step $n$.

As observed by <cit.>, bounded updates are common in practice.

§ TAIL-BOUND ANALYSIS

One application of our central-moment analysis is to use the bounds on higher (central) moments with concentration-of-measure inequalities to bound the probability that the accumulated cost deviates from some given quantity.

In this section, we sketch how we produce the tail bounds shown in <ref>(c).

There are a lot of concentration-of-measure inequalities in probability theory <cit.>.
Among those, one of the most important is Markov's inequality:

If $X$ is a nonnegative random variable and $a > 0$, then $\prob[X \ge a] \le \frac{\expe[X^k]}{a^k}$ for any $k \in \bbN$.

Recall that <ref>(b) presents upper bounds on the raw moments $\expe[\id{tick}] \le 2d+4$ and $\expe[\id{tick}^2] \le 4d^2+22d+28$ for the cost accumulator $\id{tick}$.
With Markov's inequality, we derive the following tail bounds:

\begin{align}
\prob[\id{tick} \ge 4d] & \le \frac{\expe[\id{tick}]}{4d} \le \frac{2d+4}{4d} \xrightarrow{d \to \infty} \frac{1}{2}, \label{Eq:FstMomMarkov} \\
\prob[\id{tick} \ge 4d] & \le \frac{\expe[\id{tick}^2]}{(4d)^2} \le \frac{4d^2+22d+28}{16d^2} \xrightarrow{d \to \infty} \frac{1}{4}. \label{Eq:SndMomMarkov}
\end{align}

Note that (<ref>) provides an asymptotically more precise bound on $\prob[\id{tick} \ge 4d]$ than (<ref>) does, when $d$ approaches infinity.

Central-moment analysis can obtain an even more precise tail bound.

As presented in <ref>(b), our analysis derives $\vari[\id{tick}] \le 22d+28$ for the variance of 𝕀tick.
We can now employ concentration inequalities that involve variances of random variables.
Recall Cantelli's inequality:

If $X$ is a random variable and $a > 0$, then $\prob[X \!-\! \expe[X] \!\ge\! a] \le \frac{\vari[X]}{\vari[X]+a^2}$ and $\prob[X \!-\! \expe[X] \!\le\! {-a}] \le \frac{\vari[X]}{\vari[X]+a^2}$.

With Cantelli's inequality, we obtain the following tail bound, where we assume $d \ge 2$:

\begin{equation}\label{Eq:VariCantelli}\small
\begin{split}
& \prob[\id{tick} \ge 4d] = \prob[\id{tick} - (2d+4) \ge 2d - 4] \\
{}\le{} & \prob[\id{tick} - \expe[\id{tick}] \ge 2d-4] \le \frac{\vari[\id{tick}]}{\vari[\id{tick}] + (2d-4)^2} \\
={} & 1 - \frac{(2d-4)^2}{\vari[\id{tick}] + (2d-4)^2} \le \frac{22d+28}{4d^2+ 6d+44} \xrightarrow{d \to \infty} 0.
\end{split}
\end{equation}

For all $d \ge 15$, (<ref>) gives a more precise bound than both (<ref>) and (<ref>).
It is also clear from <ref>(c), where we plot the three tail bounds (<ref>), (<ref>), and (<ref>), that the asymptotically most precise bound is the one obtained via variances.

In general, for higher central moments, we employ Chebyshev's inequality to derive tail bounds:

If $X$ is a random variable and $a>0$, then $\prob[\abs{X - \expe[X]} \ge a] \le \frac{\expe[(X-\expe[X])^{2k}]}{ a^{2k} }$ for any $k \in \bbN$.

In our experiments, we use Chebyshev's inequality to derive tail bounds from the fourth central moments.
We will show in <ref> that these tail bounds can be more precise than those obtained from both raw moments and variances.

§ IMPLEMENTATION AND EXPERIMENTS

Our tool is implemented in OCaml, and consists of about 5,300 LOC.

The tool works on imperative arithmetic probabilistic programs using a CFG-based IR <cit.>.

The language supports recursive functions, continuous distributions, unstructured control-flow, and local variables.

To infer the bounds on the central moments for a cost accumulator in a program, the user needs to specify the order of the analyzed moment, and a maximal degree for the polynomials to be used in potential-function templates.

Using APRON <cit.>, we implemented an interprocedural numeric analysis to infer the logical contexts used in the derivation.

We use the off-the-shelf solver Gurobi <cit.> for LP solving.

Evaluation setup

We evaluated our tool to answer the following three research questions:

*  How does the raw-moment inference part of our tool compare to existing techniques for expected-cost bound analysis <cit.>?

*  How does our tool compare to the state of the art in tail-probability analysis (which is based only on higher raw moments <cit.>)?

*  How scalable is our tool? Can it analyze programs with many recursive functions?

For the first question, we collected a broad suite of challenging examples from related work <cit.> with different loop and recursion patterns, as well as probabilistic branching, discrete sampling, and continuous sampling.

Our tool achieved comparable precision and efficiency with the prior work on expected-cost bound analysis <cit.>.

The details are included in
the technical report <cit.>.

For the second question, we evaluated our tool on the complete benchmarks from Kura el al. <cit.>.

We also conducted a case study of a timing-attack analysis for a program provided by DARPA during engagements of the STAC program <cit.>, where central moments are more useful than raw moments to bound the success probability of an attacker.

We include the case study in
the technical report <cit.>.

For the third question, we conducted case studies on two sets of synthetic benchmark programs:

* coupon-collector programs with $N$ coupons ($N \in [1,10^3]$), where each program is implemented as a set of tail-recursive functions, each of which represents a state of coupon collection, i.e., the number of coupons collected so far; and

* random-walk programs with $N$ consecutive one-dimensional random walks ($N \in [1,10^3]$), each of which starts at a position that equals the number of steps taken by the previous random walk to reach the ending position (the origin).

Each program is implemented as a set of non-tail-recursive functions, each of which represents a random walk.
The random walks in the same program can have different transition probabilities.

The largest synthetic program has nearly 16,000 LOC.

We then ran our tool to derive an upper bound on the fourth (resp., second) central moment of the runtime for each coupon-collector (resp., random-walk) program.

The experiments were performed on a machine with an Intel Core i7 3.6GHz processor and 16GB of RAM under macOS Catalina 10.15.7.

Some of the evaluation results to answer the second research question are presented in <ref>.

The programs (1-1) and (1-2) are coupon-collector problems with a total of two and four coupons, respectively.

The programs (2-1) and (2-2) are one-dimensional random walks with integer-valued and real-valued coordinates, respectively.

We omit three other programs here but include the full results in
the technical report <cit.>.

The table contains the inferred upper bounds on the moments for runtimes of these programs, and the running times of the analyses.
We compared our results with Kura et al.'s inference tool for raw moments <cit.>.

Our tool is as precise as, and sometimes more precise than the prior work on all the benchmark programs.
Meanwhile, our tool is able to infer an upper bound on the raw moments of degree up to four on all the benchmarks, while the prior work reports failure on some higher moments for the random-walk programs.

In terms of efficiency, our tool completed each example in less than 10 seconds, while the prior work took more than a few minutes on some programs.

One reason why our tool is more efficient is that we always reduce higher-moment inference with non-linear polynomial templates to efficient LP solving, but the prior work requires semidefinite programming (SDP) for polynomial templates.

Upper bounds on the raw/central moments of runtimes, with comparison to <cit.>. “T/O” stands for timeout after 30 minutes. “N/A” means that the tool is not applicable. “-” indicates that the tool fails to infer a bound. Entries with more precise results or less analysis time are marked in bold. Full results are included in
the technical report <cit.>.

2*program     2*moment     2c@||@this work     2c<cit.>
upper bnd.     time (sec)     upper bnd.     time (sec)
5*(1-1)     2nd raw     201     0.019     201
3rd raw     3829     0.019     3829     0.020
4th raw     90705     0.023     90705     0.027
2nd central     32     0.029     N/A     N/A
4th central     9728     0.058     N/A     N/A
5*(1-2)     2nd raw     2357     1.068     3124     0.037
3rd raw     148847     1.512     171932     0.062
4th raw      11285725     1.914     12049876     0.096
2nd central     362     3.346     N/A     N/A
4th central     955973     9.801     N/A     N/A
5*(2-1)     2nd raw     2320     0.016     2320     11.380
3rd raw     691520     0.018     -     16.056
4th raw     340107520     0.021     -     23.414
2nd central     1920     0.026     N/A     N/A
4th central     289873920     0.049     N/A     N/A
5*(2-2)     2nd raw     8375     0.022     8375     38.463
3rd raw     1362813     0.028     -     73.408
4th raw     306105209     0.035     -     141.072
2nd central      5875     0.029     N/A     N/A
4th central     447053126     0.086     N/A     N/A

Besides raw moments, our tool is also capable of inferring upper bounds on the central moments of runtimes for the benchmarks.

To evaluate the quality of the inferred central moments, <ref> plots the upper bounds of tail probabilities on runtimes $T$ obtained by <cit.>, and those by our central-moment analysis.
Specifically, the prior work uses Markov's inequality (<ref>), while we are also able to apply Cantelli's and Chebyshev's inequality (<ref>) with central moments.

Our tool outperforms the prior work on programs (1-1) and (1-2), (2-3), and (2-5), and derives better tail bounds when $d$ is large on program (2-2) and (2-4), while it obtains similar curves on program (2-1).

In <ref>, we demonstrate the running times of our tool on the two sets of synthetic benchmark programs;
<ref> plots the analysis times for coupon-collector programs as a function of the independent variable $N$ (the total number of coupons),
and <ref> plots the analysis times for random-walk programs as a function of $N$ (the total number of random walks).

The evaluation statistics show that our tool achieves good scalability in both case studies:
the runtime is almost a linear function of the program size, which is the number of recursive functions for both case studies.

Two reasons why our tool is scalable on the two sets of programs are
(i) our analysis is compositional and uses function summaries to analyze function calls, and
(ii) for a fixed set of templates and a fixed diameter of the call graph, the number of linear constraints generated by our tool grows linearly with the size of the program, and the LP solvers available nowadays can handle large problem instances efficiently.

Coupon Collector

Random Walk

Running times of our tool on two sets of synthetic benchmark programs. Each figure plots the runtimes as a function of the size of the analyzed program.

Higher central moments can also provide more information about the shape of a distribution, e.g.,

the skewness (i.e., $\frac{ \expe[(T-\expe[T])^3] }{ (\vari[T])^{\sfrac{3}{2}} }$) indicates how lopsided the distribution of $T$ is,

and the kurtosis (i.e., $\frac{ \expe[(T - \expe[T])^4 ] }{ (\vari[T])^2 }$) measures the heaviness of the tails of the distribution of $T$.

We used our tool to analyze two variants of the random-walk program (2-1).

The two random walks have different transition probabilities and step length, but they have the same expected runtime $\expe[T]$.

<ref> presents the skewness and kurtosis derived from the moment bounds inferred by our tool.

Skewness & kurtosis.

program     skewness     kurtosis
2.1362     10.5633
2.9635      17.5823

A positive skew indicates that the mass of the distribution is concentrated on the left, and larger skew means the concentration is more left.

A larger kurtosis, on the other hand, indicates that the distribution has fatter tails.

Therefore, as the derived skewness and kurtosis indicate, the distribution of the runtime $T$ for  should be more left-leaning and have fatter tails than the distribution for .

Density estimates for the runtime $T$, obtained by simulation, are shown in <ref>.

Density estimation for the runtime $T$ of two variants  and  of (2-1).

Our tool can also derive symbolic bounds on higher moments.

The table below presents the inferred upper bounds on the variances for the random-walk benchmarks, where we replace the concrete inputs with symbolic pre-conditions.

program     pre-condition     upper bnd. on the variance
(2-1)     $x \ge 0$     $1920x$
(2-2)     $x \ge 0$     $2166.6667x+1541.6667$

§ CONCLUSION

We have presented an automatic central-moment analysis for probabilistic programs that support general recursion and continuous distributions,
by deriving symbolic interval bounds on higher raw moments for cost accumulators.

We have proposed moment semirings for compositional reasoning and moment-polymorphic recursion for interprocedural analysis.

We have proved soundness of our technique using a recent extension to the Optional Stopping Theorem.

The effectiveness of our technique has been demonstrated with our prototype implementation and the analysis of a broad suite of benchmarks, as well as a tail-bound analysis.

In the future, we plan to go beyond arithmetic programs and add support for more datatypes, e.g., Booleans and lists.

We will also work on other kinds of uncertain quantities for probabilistic programs.

Another research direction is to apply our analysis to higher-order functional programs.

This article is based on research supported, in part, by a gift from Rajiv and Ritu Batra;
by ONR under grants N00014-17-1-2889 and N00014-19-1-2318;
by DARPA under AA contract FA8750-18-C0092;
and by the NSF under SaTC award 1801369, SHF awards 1812876 and 2007784, and CAREER
award 1845514.
Any opinions, findings, and conclusions or recommendations
expressed in this publication are those of the authors,
and do not necessarily reflect the views of the sponsoring

We thank Jason Breck and John Cyphert for formulating, and providing
us with an initial formalization of the problem of timing-attack

§ PRELIMINARIES ON MEASURE THEORY

Interested readers can refer to textbooks and notes in the literature <cit.> for more details.

§.§ Basics

A measurable space is a pair $(S,\calS)$, where $S$ is a nonempty set, and $\calS$ is a $\sigma$-algebra on $S$, i.e., a family of subsets of $S$ that contains $\emptyset$ and is closed under complement and countable unions.

The smallest $\sigma$-algebra that contains a family $\calA$ of subsets of $S$ is said to be generated by $\calA$, denoted by $\sigma(\calA)$.

Every topological space $(S,\tau)$ admits a Borel $\sigma$-algebra, given by $\sigma(\tau)$.

This gives canonical $\sigma$-algebras on $\bbR$, $\bbQ$, $\bbN$, etc.

A function $f : S \to T$, where $(S,\calS)$ and $(T,\calT)$ are measurable spaces, is said to be $(\calS,\calT)$-measurable, if $f^{-1}(B) \in \calS$ for each $B \in \calT$.

If $T = \bbR$, we tacitly assume that the Borel $\sigma$-algebra is defined on $T$, and we simply call $f$ measurable, or a random variable.

Measurable functions form a vector space, and products, maxima, and limiting operations preserve measurability.

A measure $\mu$ on a measurable space $(S,\calS)$ is a mapping from $\calS$ to $[0,\infty]$ such that
(i) $\mu(\emptyset) = 0$, and
(ii) for all pairwise-disjoint $\{A_n\}_{n \in \bbZ^+}$ in $\calS$, it holds that $\mu(\bigcup_{n \in \bbZ^+} A_i) = \sum_{n \in \bbZ^+} \mu(A_i)$.

The triple $(S,\calS,\mu)$ is called a measure space.

A measure $\mu$ is called a probability measure, if $\mu(S) = 1$.

We denote the collection of probability measures on $(S,\calS)$ by $\bbD(S,\calS)$.

For each $x \in S$, the Dirac measure $\delta(x)$ is defined as $\lambda A. [x \in A]$.

For measures $\mu$ and $\nu$, we write $\mu + \nu$ for the measure $\lambda A. \mu(A) + \nu(A)$.

For measure $\mu$ and scalar $c \ge 0$, we write $c \cdot \mu$ for the measure $\lambda A. c \cdot \mu(A)$.

The integral of a measurable function $f$ on $A \in \calS$ with respect to a measure $\mu$ on $(S,\calS)$ is defined following Lebesgue's theory and is denoted by $\mu(f;A)$, $\int_A f d\mu$, or $\int_A f(x) \mu(dx)$.

If $\mu$ is a probability measure, we call the integral as the expectation of $f$, written $\expe_{x \sim \mu}[f; A]$, or simply $\expe[f; A]$ when the scope is clear in the context.

If $A = S$, we tacitly omit $A$ from the notations.

For each $A \in \calS$, it holds that $\mu(f; A) = \mu(f \mathrm{I}_A)$, where $\mathrm{I}_A$ is the indicator function for $A$.

We denote the collection of integrable functions on $(S,\calS,\mu)$ by $\calL^1(S,\calS,\mu)$.

A kernel from a measurable space $(S,\calS)$ to another $(T,\calT)$ is a mapping from $S \times \calT$ to $[0,\infty]$ such that:
(i) for each $x \in S$, the function $\lambda B. \kappa(x,B)$ is a measure on $(T,\calT)$, and
(ii) for each $B \in \calT$, the function $\lambda x. \kappa(x,B)$ is measurable.

We write $\kappa : (S,\calS) \rightsquigarrow (T,\calT)$ to declare that $\kappa$ is a kernel from $(S,\calS)$ to $(T,\calT)$.

Intuitively, kernels describe measure transformers from one measurable space to another.

A kernel $\kappa$ is called a probability kernel, if $\kappa(x,T) = 1$ for all $x \in S$.

We denote the collection of probability kernels from $(S,\calS)$ to $(T,\calT)$ by $\bbK((S,\calS),(T,\calT))$.

If the two measurable spaces coincide, we simply write $\bbK(S,\calS)$.

We can “push-forward” a measure $\mu$ on $(S,\calS)$ to a measure on $(T,\calT)$ through a kernel $\kappa: (S,\calS) \rightsquigarrow (T,\calT)$ by integration:[
We use a monad bind notation $\bind$ here.

Indeed, the category of measurable spaces admits a monad with sub-probability measures <cit.>.
\mu \bind \kappa \defeq  \lambda B. \int_S  \kappa(x,B) \mu(dx).

§.§ Product Measures

The product of two measurable spaces $(S,\calS)$ and $(T,\calT)$ is defined as $(S,\calS) \otimes (T,\calT) \defeq (S \times T, \calS \otimes \calT)$, where $\calS \otimes \calT$ is the smallest $\sigma$-algebra that makes coordinate maps measurable, i.e., $\sigma( \{ \pi_1^{-1}(A) \mid A \in \calS\} \cup \{ \pi_2^{-1}(B) \mid B \in \calT \}  ) )$, where $\pi_i$ is the $i$-the coordinate map.

If $\mu_1$ and $\mu_2$ are two probability measures on $(S,\calS)$ and $(T,\calT)$, respectively, then there exists a unique probability measure $\mu$ on $(S,\calS) \otimes (T,\calT)$, by Fubini's theorem, called the product measure of $\mu_1$ and $\mu_2$, written $\mu_1 \otimes \mu_2$, such that $\mu(A \times B) = \mu_1(A) \cdot \mu_2(B)$ for all $A \in \calS$ and $B \in \calT$.

If $\mu$ is a probability measure on $(S,\calS)$ and $\kappa : (S,\calS) \rightsquigarrow (T,\calT)$ is a probability kernel, then we can construct the a probability measure on $(S,\calS) \otimes (T,\calT)$ that captures all transitions from $\mu$ through $\kappa$:
\mu \otimes \kappa \defeq \lambda(A,B). \int_A \kappa(x,B) \mu (dx).

If $\mu$ is a probability measure on $(S_0,\calS_0)$ and $\kappa_i : (S_{i-1},\calS_{i-1}) \rightsquigarrow (S_i,\calS_i)$ is a probability kernel for $i=1,\cdots,n$,where $n \in \bbZ^+$, then we can construct a probability measure on $\bigotimes_{i=0}^n (S_i,\calS_i)$, i.e., the space of sequences of $n$ transitions by iteratively applying the kernels to $\mu$:

\begin{align*}
\mu \otimes \bigotimes\nolimits_{i=1}^{0} \kappa_i & \defeq \mu, \\
\mu \otimes \bigotimes\nolimits_{i=1}^{k+1} \kappa_i & \defeq (\mu \otimes \bigotimes\nolimits_{i=1}^{k} \kappa_i) \otimes \kappa_{k+1}, & 0 \le k < n.
\end{align*}

Let $(S_i,\calS_i)$, $i \in \calI$ be a family of measurable spaces.

Their product, denoted by $\bigotimes_{i \in \calI} (S_i,\calS_i) = (\prod_{i \in \calI} S_i, \bigotimes_{i \in \calI} \calS_i)$, is the product space with the smallest $\sigma$-algebra such that for each $i \in \calI$, the coordinate map $\pi_i$ is measurable.

The theorem below is widely used to construct a probability measures on an infinite product via probability kernels.

Let $(S_i,\calS_i)$, $i \in \bbZ^+$ be a sequence of measurable spaces.

Let $\mu_0$ be a probability measure on $(S_0,\calS_0)$.

For each $i \in \bbN$, let $\kappa_i : \bigotimes_{k=0}^{i-1} (S_k,\calS_k) \rightsquigarrow (S_i,\calS_i)$ be a probability kernel.

Then there exists a sequence of probability measures $\mu_i \defeq \mu_0 \otimes \bigotimes_{k=1}^i \kappa_k$, $i \in \bbZ^+$, and there exists a uniquely defined probability measure $\mu$ on $\bigotimes_{k=0}^\infty (S_k,\calS_k)$ such that $\mu_i(A) = \mu(A \times \prod_{k=i+1}^\infty S_k)$ for all $i \in \bbZ^+$ and $A \in \bigotimes_{k=0}^i \calS_i$.

§.§ Conditional Expectations

Let $(\Omega,\calF,\prob)$ be a measure space, where $\prob$ is a probability measure.

We write $\calL^1$ for $\calL^1(\Omega,\calF,\prob)$.

Let $X \in \calL^1$ be an integrable random variable, and $\calG \subseteq \calF$ be a sub-$\sigma$-algebra of $\calF$.

Then there exists a random variable $Y \in \calL^1$ such that:
(i) $Y$ is $\calG$-measurable, and
(ii) for each $G \in \calG$, it holds that $\expe[Y; G] = \expe[X; G]$.

Such a random variable $Y$ is said to be a version of conditional expectation $\expe[X \mid \calG]$ of $X$ given $\calG$.

Conditional expectations admit almost-sure uniqueness.

On the other hand, if $X$ is a nonnegative random variable, the conditional expectation $\expe[X \mid \calG]$ is also defined and almost-surely unique.

Intuitively, for some $\omega \in \Omega$, $\expe[X \mid \calG](\omega)$ is the expectation of $X$ given the set of values $Z(\omega)$ for every $\calG$-measurable random variable $Z$.

For example, if $\calG = \{\emptyset,\Omega\}$, which contains no information, then $\expe[X \mid \calG](\omega) = \expe[X]$.

If $\calG = \calF$, which contains full information, then $\expe[X \mid \calG](\omega) = X(\omega)$.

We review some useful properties of conditional expectations below.

Let $X \in \calL^1$, and $\calG,\calH$ be sub-$\sigma$-algebras of $\calF$.

* If $Y$ is any version of $\expe[X \mid \calG]$, then $\expe[Y] = \expe[X]$.

* If $X$ is $\calG$-measurable, then $\expe[X \mid \calG] = X$, a.s.

* If $\calH$ is a sub-$\sigma$-algebra of $\calG$, then $\expe[\expe(X \mid \calG) \mid \calH] = \expe[X \mid \calH]$, a.s.

* If $Z$ is $\calG$-measurable and $Z \cdot X \in \calL^1$, then $\expe[Z \cdot X \mid \calG] = Z \cdot \expe[X \mid \calG]$, a.s.

On the other hand, if $X$ is a nonnegative random variable and $Z$ is a nonnegative $\calG$-measurable random variable, then the property also holds, with the understanding that both sides might be infinite.

§.§ Convergence Theorems

We review two important convergence theorems for series of random variables.

If $\{f_n\}_{n \in \bbZ^+}$ is a non-decreasing sequence of nonnegative measurable functions on a measure space $(S,\calS,\mu)$, and $\{f_n\}_{n \in \bbZ^+}$ converges to $f$ pointwise, then $f$ is measurable and
\lim_{n \to \infty} \mu(f_n) = \mu(f) \le \infty.

Further, the theorem still holds if $f$ is chosen as a measurable function and “$\{f_n\}_{n \in \bbZ^+}$ converges to $f$ pointwise” holds almost everywhere, rather than everywhere.

If $\{f_n\}_{n \in \bbZ^+}$ is a sequence of measurable functions on a measure space $(S,\calS,\mu)$, $\{f_n\}_{n \in \bbZ^+}$ converges to $f$ pointwise, and $\{f_n\}_{n \in \bbZ^+}$ is dominated by a nonnegative integrable function $g$ (i.e., $|f_n(x)| \le g(x)$ for all $n \in \bbZ^+$, $x \in S$), then $f$ is integrable and
\lim_{n \to \infty} \mu(f_n) = \mu(f).

Further, the theorem still holds if $f$ is chosen as a measurable function and “$\{f_n\}_{n \in \bbZ^+}$ converges to $f$ pointwise and is dominated by $g$” holds almost everywhere.

§ OPERATIONAL COST SEMANTICS

We follow a distribution-based approach <cit.> to define reduction rules for the semantics of .

A probabilistic semantics steps a program configuration to a probability distribution on configurations.

To formally describe these distributions, we need to construct a measurable space of program configurations.

Our approach is to construct a measurable space for each of the four components of configurations, and then use their product measurable space as the semantic domain.

* Valuations $\gamma : \mathsf{VID} \to \bbR$ are finite real-valued maps, so we define $(V,\calV) \defeq (\bbR^{\mathsf{VID}},\calB(\bbR^{\mathsf{VID}}))$ as the canonical structure on a finite-dimensional space.

* The executing statement $S$ can contain real numbers, so we need to “lift” the Borel $\sigma$-algebra on $\bbR$ to program statements.

Intuitively, statements with exactly the same structure can be treated as vectors of parameters that correspond to their real-valued components.

Formally, we achieve this by constructing a metric space on statements and then extracting a Borel $\sigma$-algebra from the metric space.

<ref> presents a recursively defined metric $d_S$ on statements, as well as metrics $d_E$, $d_L$, and $d_D$ on expressions, conditions, and distributions, respectively, as they are required by $d_S$.

We denote the result measurable space by $(S,\calS)$.

* A continuation $K$ is either an empty continuation $\kstop$, a loop continuation $\kloop{L}{S}{K}$, or a sequence continuation $\kseq{S}{K}$.

Similarly, we construct a measurable space $(K,\calK)$ on continuations by extracting from a metric space.

<ref> shows the definition of a metric $d_K$ on continuations.

* The cost accumulator $\alpha \in \bbR$ is a real number, so we define $(W,\calW) \defeq (\bbR,\calB(\bbR))$ as the canonical measurable space on $\bbR$.

Then the semantic domain is defined as the product measurable space of the four components: $(\Sigma,\calO) \defeq (V,\calV) \otimes (S,\calS) \otimes (K,\calK) \otimes (W,\calW)$.

\begin{align*}
d_E(x,x) & \defeq 0 \\
d_E(c_1,c_2) & \defeq \abs{c_1-c_2} \\
d_E(\eadd{E_{11}}{E_{12}},\eadd{E_{21}}{E_{22}}) & \defeq d_E(E_{11},E_{21}) + d_E(E_{12},E_{22}) \\
d_E(\emul{E_{11}}{E_{12}},\emul{E_{21}}{E_{22}}) & \defeq d_E(E_{11},E_{21}) + d_E(E_{12},E_{22}) \\
d_E(E_1,E_2) & \defeq \infty~\text{otherwise}
\end{align*}

\begin{align*}
d_L(\ctop,\ctop) & \defeq 0 \\
d_L(\cneg{L_1},\cneg{L_2}) & \defeq d_L(L_1,L_2) \\
d_L(\cconj{L_{11}}{L_{12}},\cconj{L_{21}}{L_{22}}) & \defeq d_L(L_{11},L_{21}) + d_L(L_{12},L_{22}) \\
d_L(\cle{E_{11}}{E_{12}},\cle{E_{21}}{E_{22}}) & \defeq d_E(E_{11},E_{21}) + d_E(E_{12},E_{22}) \\
d_L(L_1,L_2) & \defeq \infty~\text{otherwise}
\end{align*}

\begin{align*}
d_D(\kw{uniform}(a_1,b_1),\kw{uniform}(a_2,b_2)) & \defeq |a_1-a_2| + |b_1-b_2| \\
d_D(D_1,D_2) & \defeq \infty~\text{otherwise}
\end{align*}

\begin{align*}
d_S(\iskip,\iskip) & \defeq 0 \\
d_S(\itick{c_1},\itick{c_2}) & \defeq |c_1 - c_2| \\
d_S(\iassign{x}{E_1},\iassign{x}{E_2}) & \defeq d_E(E_1,E_2) \\
d_S(\isample{x}{D_1},\isample{x}{D_2}) & \defeq d_D(D_1,D_2) \\
d_S(\iinvoke{f},\iinvoke{f} ) & \defeq 0 \\
d_S(\iprob{p_1}{S_{11}}{S_{12}},\iprob{p_2}{S_{21}}{S_{22}}) & \defeq |p_1-p_2| + d_S(S_{11},S_{21}) + d_S(S_{12},S_{22}) \\
d_S(\icond{L_1}{S_{11}}{S_{12}},\icond{L_2}{S_{21}}{S_{22}}) & \defeq d_L(L_1,L_2) + d_S(S_{11},S_{21}) + d_S(S_{12},S_{22}) \\
d_S(\iloop{L_1}{S_1}, \iloop{L_2}{S_2}) & \defeq d_L(L_1,L_2)+d_S(S_1,S_2) \\
d_S(S_{11};S_{12} , S_{21};S_{22}) & \defeq d_S(S_{11},S_{21}) + d_S(S_{12},S_{22}) \\
d_S(S_1,S_2) & \defeq \infty ~\text{otherwise}
\end{align*}

\begin{align*}
d_K(\kstop, \kstop ) & \defeq 0 \\
d_K(\kloop{L_1}{S_1}{K_1},\kloop{L_2}{S_2}{K_2}) & \defeq d_L(L_1,L_2) + d_S(S_1,S_2) + d_K(K_1,K_2) \\
d_K(\kseq{S_1}{K_1},\kseq{S_2}{K_2}) & \defeq d_S(S_1,S_2) + d_K(K_1,K_2) \\
d_K(K_1,K_2) & \defeq \infty ~\text{otherwise}
\end{align*}

Metrics for expressions, conditions, distributions, statements, and continuations.

<ref> presents the rules of the evaluation relation $\sigma \mapsto \mu$ for  where $\sigma$ is a configuration and $\mu$ is a probability distribution on configurations.

Note that in , expressions $E$ and conditions $L$ are deterministic, so we define a standard big-step evaluation relation for them, written $\gamma \vdash E \Downarrow r$ and $\gamma \vdash L \Downarrow b$, where $\gamma$ is a valuation, $r \in \bbR$, and $b \in \{ \top, \bot \}$.

Most of the rules in <ref>, except (E-Sample) and (E-Prob), are also deterministic as they step to a Dirac measure.

$\gamma \vdash E \Downarrow r $    “the expression $E$ evaluates to a real value $r$ under the valuation $\gamma$”
E-Var γ(x) = r  γ⊢x ⇓r
E-Const  γ⊢c ⇓c
E-Add γ⊢E_1 ⇓r_1
γ⊢E_2 ⇓r_2
r = r_1 + r_2  γ⊢ E_1  E_2 ⇓r
E-Mul γ⊢E_1 ⇓r_1
γ⊢E_2 ⇓r_2
r = r_1 ·r_2γ⊢E_1 E_2 ⇓r

$\gamma \vdash L \Downarrow b$    “the condition $L$ evaluates to a Boolean value $b$ under the valuation $\gamma$”
E-Top  γ⊢⇓⊤
E-Neg γ⊢L ⇓b  γ⊢L ⇓b
E-Conj γ⊢L_1 ⇓b_1
γ⊢L_2 ⇓b_2  γ⊢L_1L_2 ⇓b_1 ∧b_2
E-Le γ⊢E_1 ⇓r_1
γ⊢E_2 ⇓r_2  γ⊢E_1E_2 ⇓[r_1 ≤r_2]

$\tuple{\gamma,S,K,\alpha} \mapsto \mu$    “the configuration $\tuple{\gamma,S,K,\alpha}$ steps to a probability distribution $\mu$ on $\tuple{\gamma',S',K',\alpha'}$'s”
E-Skip-Stop  γ,,,α ↦δ(γ,,,α)
E-Skip-Loop γ⊢L ⇓b  γ, , SLK ,α ↦[b] ·δ(γ,S, SLK ,α) + [b] ·δ( γ,, K,α )
E-Skip-Seq  γ,,SK,α ↦δ( γ,S,K,α )
E-Tick   γ,c  ,K,α ↦δ(γ, ,K,α+c )
E-Assign γ⊢E ⇓r  γ, xE ,K,α ↦δ( γ[x ↦r], ,K,α )
E-Sample   γ,xD ,K,α ↦μ_D λr. δ(  γ[x ↦r],  ,K,α )
E-Call   γ, f  ,K,α ↦δ( γ, (f) , K,α )
E-Prob  γ, pS_1S_2 ,K,α ↦p ·δ(γ,S_1,K,α) + (1-p) ·δ(γ,S_2,K,α)
E-Cond  γ⊢L ⇓b  γ,LS_1S_2 ,K,α ↦[b] ·δ(γ,S_1,K,α) + [b] ·δ(γ,S_2,K,α)
E-Loop  γ,LS,K,α ↦δ(  γ,,LSK,α )
E-Seq   γ,S_1;S_2, K,α ↦δ(  γ,S_1,S_2K,α )

Rules of the operational semantics of .

The evaluation relation $\mapsto$ can be interpreted as a distribution transformer.

Indeed, $\mapsto$ can be seen as a probability kernel.

Let $\gamma :\mathsf{VID} \to \bbR$ be a valuation.

* Let $E$ be an expression. Then there exists a unique $r \in \bbR$ such that $\gamma \vdash E \Downarrow r$.

* Let $L$ be a condition. Then there exists a unique $b \in \{ \top, \bot \}$ such that $\gamma \vdash L \Downarrow b$.

By induction on the structure of $E$ and $L$.

For every configuration $\sigma \in \Sigma$, there exists a unique $\mu\in\bbD(\Sigma,\calO)$ such that $\sigma \mapsto \mu$.

Let $\sigma = \tuple{\gamma,S,K,\alpha}$.
Then by case analysis on the structure of $S$, followed by a case analysis on the structure of $K$ if $S = \iskip$.
The rest of the proof appeals to <ref>.

The evaluation relation $\mapsto$ defines a probability kernel on program configurations.

<ref> tells us that $\mapsto$ can be seen as a function $\hat{\mapsto}$ defined as follows:

\[
\hat{\mapsto}(\sigma,A) \defeq \mu(A) \quad \text{where} \quad \sigma \mapsto \mu.
\]

It is clear that $\lambda A. \hat{\mapsto}(\sigma,A)$ is a probability measure.

On the other hand, to show that $\lambda\sigma. \hat{\mapsto}(\sigma,A)$ is measurable for any $A \in \calO$, we need to prove that for $B \in \calB(\bbR)$, it holds that $\scrO(A,B) \defeq (\lambda\sigma. \hat{\mapsto}(\sigma,A))^{-1}(B) \in \calO$.

We introduce skeletons of programs to separate real numbers and discrete structures.

\begin{align*}
\hat{S} & \Coloneqq \iskip \mid \itick{\square_\ell} \mid \iassign{x}{\hat{E}} \mid \isample{x}{\hat{D}} \mid \iinvoke{f} \\
& \mid \iprob{\square_\ell}{\hat{S}_1}{\hat{S}_2} \mid \icond{\hat{L}}{\hat{S}_1}{\hat{S}_2} \\
& \mid \iloop{\hat{L}}{\hat{S}} \mid \hat{S}_1; \hat{S}_2 \\
\hat{L} & \Coloneqq \ctop \mid \cneg{\hat{L}} \mid \cconj{\hat{L}_1}{\hat{L}_2} \mid \cle{\hat{E}_1}{\hat{E}_2} \\
\hat{E} & \Coloneqq x \mid \square_\ell \mid \eadd{\hat{E}_1}{\hat{E}_2} \mid \emul{\hat{E}_1}{\hat{E}_2} \\
\hat{D} & \Coloneqq \kw{uniform}(\square_{\ell_a},\square_{\ell_b}) \\
\hat{K} & \Coloneqq \kstop \mid \kloop{\hat{L}}{\hat{S}}{\hat{K}} \mid \kseq{\hat{S}}{\hat{K}}
\end{align*}

The holes $\square_\ell$ are placeholders for real numbers parameterized by locations $\ell \in \m{LOC}$.

We assume that the holes in a program structure are always pairwise distinct.

Let $\eta : \m{LOC} \to \bbR$ be a map from holes to real numbers and $\eta(\hat{S})$ (resp., $\eta(\hat{L})$, $\eta(\hat{E})$, $\eta(\hat{D})$, $\eta(\hat{K})$) be the instantiation of a statement (resp., condition, expression, distribution, continuation) skeleton by substituting $\eta(\ell)$ for $\square_\ell$.

One important property of skeletons is that the “distance” between any concretizations of two different skeletons is always infinity with respect to the metrics in <ref>.

Observe that

\[
\scrO(A,B) = \bigcup_{\hat{S},\hat{K}} \scrO(A,B) \cap \{ \tuple{\gamma,\eta(\hat{S}),\eta(\hat{K}),\alpha} \mid \text{any}~\gamma,\alpha,\eta  \}
\]

and that $\hat{S},\hat{K}$ are countable families of statement and continuation skeletons.

Thus it suffices to prove that every set in the union, which we denote by $\scrO(A,B) \cap \scrC(\hat{S},\hat{K})$ later in the proof, is measurable.

Note that $\scrC(\hat{S},\hat{K})$ itself is indeed measurable.

Further, the skeletons $\hat{S}$ and $\hat{K}$ are able to determine the evaluation rule for all concretized configurations.

Thus we can proceed by a case analysis on the evaluation rules.

To aid the case analysis, we define a deterministic evaluation relation $\xmapsto{\text{det}}$ by getting rid of the $\delta(\cdot)$ notations in the rules in <ref> except probabilistic ones (E-Sample) and (E-Prob).

Obviously, $\xmapsto{\text{det}}$ can be interpreted as a measurable function on configurations.

* If the evaluation rule is deterministic, then we have

\begin{align*}
& \scrO(A,B) \cap \scrC(\hat{S},\hat{K}) \\
={} & \{ \sigma \mid \sigma \mapsto \mu, \mu(A) \in B \} \cap \scrC(\hat{S},\hat{K}) \\
={} & \{ \sigma \mid \sigma \xmapsto{\text{det}} \sigma', [\sigma' \in A] \in B \} \cap \scrC(\hat{S},\hat{K}) \\
={} &
\begin{dcases*}
\scrC(\hat{S},\hat{K}) & if $\{0,1\} \subseteq B$ \\
{\xmapsto{\text{det}}}^{-1}(A) \cap \scrC(\hat{S},\hat{K}) & if $1 \in B$ and $0 \not\in B$ \\
{\xmapsto{\text{det}}}^{-1}(A^c) \cap \scrC(\hat{S},\hat{K}) & if $0 \in B$ and $1 \not\in B$ \\
\emptyset & if $\{0,1\} \cap B = \emptyset$.
\end{dcases*}
\end{align*}

The sets in all the cases are measurable, so is the set $\scrO(A,B) \cap \scrC(\hat{S},\hat{K})$.

* (E-Prob): Consider $B$ with the form $(-\infty,t]$ with $t \in \bbR$.

If $t \ge 1$, then $\scrO(A,B) = \Sigma$.

Otherwise, let us assume $t < 1$.

Let $\hat{S} = \iprob{\square}{\hat{S}_1}{\hat{S}_2}$.

Then we have

\begin{align*}
& \scrO(A,B) \cap \scrC(\hat{S},\hat{K}) \\
={} & \{ \sigma \mid \sigma \mapsto \mu, \mu(A) \in B \} \cap \scrC(\hat{S},\hat{K}) \\
={} & \{ \sigma \mid \sigma \mapsto p \cdot \delta(\sigma_1) + (1-p) \cdot \delta(\sigma_2),\\
& \quad p \cdot [\sigma_1 \in A] + (1-p) \cdot [\sigma_2 \in A] \in B \} \cap \scrC(\hat{S},\hat{K}) \\
={} & \scrC(\hat{S},\hat{K}) \cap
\{  \tuple{\gamma,\iprob{p}{S_1}{S_2},K,\alpha} \mid \\
& \quad p \cdot [\tuple{\gamma,S_1,K,\alpha} \in A] + (1-p) \cdot [\tuple{\gamma,S_2,K,\alpha} \in A] \le t  \}  \\
={} & \scrC(\hat{S},\hat{K}) \cap {} \\
& \{ \tuple{\gamma,\iprob{p}{S_1}{S_2},K,\alpha} \mid \\
& \quad (\tuple{\gamma,S_1,K,\alpha} \in A, \tuple{\gamma,S_2,K,\alpha}\not\in A, p \le t)  \vee {} \\
& \quad (\tuple{\gamma,S_2,K,\alpha} \in A, \tuple{\gamma,S_1,K,\alpha}\not\in A, 1 - p \le t) \}.
\end{align*}

The set above is measurable because $A$ and $A^c$ are measurable, as well as $\{ p \le t\}$ and $\{ p \ge 1 - t \}$ are measurable in $\bbR$.

* (E-Sample): Consider $B$ with the form $(-\infty,t]$ with $t \in \bbR$.

Similar to the previous case, we assume that $t < 1$.

Let $\hat{S} = \isample{x}{\kw{uniform}(\square_{\ell_a},\square_{\ell_b})}$, without loss of generality.

Then we have

\begin{align*}
& \scrO(A,B) \cap \scrC(\hat{S},\hat{K}) \\
={} & \{\sigma \mid \sigma \mapsto \mu, \mu(A) \in B \} \cap \scrC(\hat{S},\hat{K}) \\
={} & \{ \sigma \mid \sigma \mapsto \mu_D \bind \kappa_{\sigma}, \int \kappa_{\sigma}(r)(A) \mu_D(dr) \le t \} \cap \scrC(\hat{S},\hat{K}) \\
={} & \scrC(\hat{S},\hat{K}) \cap  \{ \sigma \mid \sigma \mapsto \mu_D \bind \kappa_{\sigma}, \\
& \quad \mu_D(\{ r \mid \tuple{\gamma[x \mapsto r], \iskip,K,\alpha} \in A \}) \le t\} \\
={} & \scrC(\hat{S},\hat{K}) \cap \{ \tuple{\gamma,\isample{x}{\kw{uniform}(a,b)},K,\alpha} \mid a<b, \\
& \quad \mu_{\kw{uniform}(a,b)}(\{ r \mid \tuple{\gamma[x \mapsto r], \iskip,K,\alpha} \in A \}) \le t  \},
\end{align*}

where $\kappa_{\tuple{\gamma, S ,K,\alpha}} \defeq \lambda r. \delta(\tuple{\gamma[x \mapsto r], \iskip, K, \alpha})$.

For fixed $\gamma,K,\alpha$, the set $\{r \mid \tuple{\gamma[x \mapsto r],\iskip,K,\alpha} \in A\}$ is measurable in $\bbR$.

For the distributions considered in this paper, there is a sub-probability kernel $\kappa_D : \bbR^{\mathrm{ar}(D)} \rightsquigarrow \bbR$.

For example, $\kappa_{\kw{uniform}}(a,b)$ is defined to be $\mu_{\kw{uniform}(a,b)}$ if $a < b$, or $\mathbf{0}$ otherwise.

Therefore, $\lambda(a,b). \kappa_{\kw{uniform}}(a,b)(\{r \mid \tuple{\gamma[x \mapsto r],\iskip,K,\alpha} \in A\})$ is measurable, and its inversion on $(-\infty,t]$ is a measurable set on distribution parameters $(a,b)$.

Hence the set above is measurable.

§ TRACE-BASED COST SEMANTICS

To reason about moments of the accumulated cost, we follow the Markov-chain-based reasoning <cit.> to develop a trace-based cost semantics for .

Let $(\Omega,\calF) \defeq \bigotimes_{n=0}^\infty (\Sigma,\calO)$ be a measurable space of infinite traces on program configurations.

Let $\{\calF_n\}_{n \in \bbZ^+}$ be a filtration, i.e., an increasing sequence $\calF_0 \subseteq \calF_1 \subseteq \cdots \subseteq \calF$ of sub-$\sigma$-algebras in $\calF$, generated by coordinate maps $X_n(\omega) \defeq \omega_n$ for $n \in \bbZ^+$.

Let $\tuple{\scrD,S_{\mathsf{main}}}$ be an  program.

Let $\mu_0 \defeq \delta(\tuple{ \lambda\_.0 ,S_{\mathsf{main}},\kstop ,0 })$ be the initial distribution.

Let $\prob$ be the probability measure on infinite traces induced by <ref> and <ref>.

Then $(\Omega,\calF,\prob)$ forms a probability space on infinite traces of the program.

§ TRACE-BASED REASONING ON EXPECTATIONS

Recall that we define a stopping time $T : \Omega \to \bbZ^+ \cup \{\infty\}$ as

\[
T(\omega) \defeq \inf\{n \in \bbZ^+ \mid \omega_n = \tuple{\_,\iskip,\kstop,\_} \},
\]

and random variables $\{A_n\}_{n \in \bbZ^+}$, $\{\Phi_n\}_{n \in \bbZ^+}$, $\{Y_n\}_{n \in \bbZ^+}$ as

\begin{align*}
A_n(\omega) \defeq \alpha_n,
\Phi_n(\omega)  \defeq \phi(\omega_n),
Y_n(\omega) \defeq A_n(\omega) + \Phi_n(\omega),
\end{align*}

where $\omega_n=\tuple{\_,\_,\_,\alpha_n}$, and $\phi : \Sigma \to \bbR$ is an expected-potential function for expected-cost bound analysis.

Taking the stopping time into consideration, we define the stopped version for these random variables as $A_T(\omega) \defeq A_{T(\omega)}(\omega)$, $\Phi_T(\omega) \defeq \Phi_{T(\omega)}(\omega)$, $Y_T(\omega) \defeq Y_{T(\omega)}(\omega)$.

Note that $A_T,\Phi_T,Y_T$ are not guaranteed to be well-defined everywhere.

If $\prob[T < \infty] = 1$, i.e., the program terminates almost surely, then $A_T$ is well-defined almost surely and $\prob[\lim_{n \to \infty} A_n = A_T] = 1$.

Further, if $\{ A_n\}_{n \in \bbZ^+}$ is pointwise non-decreasing, then $\lim_{n \to \infty} \expe[A_n] = \expe[A_T]$.

By the property of the operational semantics, for $\omega \in \Omega$ such that $T(\omega)<\infty$, we have $A_n(\omega) = A_T(\omega)$ for all $n \ge T(\omega)$.

Then we have

\begin{align*}
& \prob[\lim_{n \to \infty} A_n = A_T] \\
={} & \prob( \{ \omega \mid \lim_{n \to \infty} A_n(\omega) = A_T(\omega) \}) \\
{}\ge{} & \prob(\{ \omega \mid \lim_{n \to \infty} A_n(\omega) = A_T(\omega) \wedge T(\omega) < \infty\}) \\
={} & \prob(\{ \omega \mid A_{T(\omega)}(\omega) = A_T(\omega) \wedge T(\omega) < \infty \}) \\
={} & \prob(\{ \omega \mid T(\omega) < \infty \}) \\
={} & 1.
\end{align*}

Now let us assume that $\{A_n\}_{n \in \bbZ^+}$ is pointwise non-decreasing.

By the property of the operational semantics, we know that $A_0 = 0$.

Therefore, $A_n$'s are nonnegative random variables, and their expectations $\expe[A_n]$'s are well-defined.

By <ref>, we have $\lim_{n \to \infty} \expe[A_n] = \expe[\lim_{n \to \infty} A_n]$.

We then conclude by the fact that $\lim_{n \to \infty} A_n = A_T$, a.s., which we just proved.

We reformulate the martingale property using the filtration $\{\calF_n\}_{n \in \bbZ^+}$.

For all $n \in \bbZ^+$, it holds that

\[
\expe[Y_{n+1} \mid \calF_n] = Y_n, \text{a.s.},
\]

i.e., the expectation of $Y_n$ conditioned on the execution history is an invariant for $n \in \bbZ^+$.

We say that a sequence of random variables $\{X_n\}_{n \in \bbZ^+}$ is adapted to a filtration $\{\calF_n\}_{n \in \bbZ^+}$ if for each $n \in \bbZ^+$, $X_n$ is $\calF_n$-measurable.

Then $\{\Phi_n\}_{n \in \bbZ^+}$ and $\{A_n\}_{n \in\bbZ^+}$ are adapted to the coordinate-generated filtration $\{\calF_n\}_{n \in \bbZ^+}$ as $\Phi_n(\omega)$ and $A_n(\omega)$ depend on $\omega_n$.

Then we have

\begin{align*}
& \expe[Y_{n+1} \mid \calF_n](\omega) \\
={} & \expe[A_{n+1} + \Phi_{n+1} \mid \calF_n](\omega) \\
={} & \expe[(A_{n+1} - A_n) +\Phi_{n+1} + A_n \mid \calF_n](\omega) \\
={} & \expe[(A_{n+1} - A_n) + \Phi_{n+1} \mid \calF_n](\omega) + A_n(\omega) \\
={} & \expe[(\alpha_{n+1} - \alpha_n) + \phi(\omega_{n+1}) \mid \calF_n] + A_n(\omega) \\
={} & \expe_{\sigma' \sim {\mapsto}(\omega_n)}[(\alpha'-\alpha_n) + \phi(\sigma')] + A_n(\omega) \\
={} & \phi(\omega_n) + A_n(\omega) \\
={} & \Phi_n(\omega) + A_n(\omega) \\
={} & Y_n(\omega), \text{a.s.}
\end{align*}

Furthermore, we have the following corollary:

\[
\expe[Y_{n+1}] = \expe[\expe(Y_{n+1} \mid \calF_n)] = \expe[Y_n],
\]

for each $n \in \bbZ^+$, thus $\expe[Y_n] = \expe[Y_0]$ for all $n \in \bbZ^+$.

Now we can prove soundness of the extended OST.

If $\expe[|Y_n|] < \infty$ for all $n \in \bbZ^+$, then $\expe[Y_T]$ exists and $\expe[Y_T] = \expe[Y_0]$ in the following situation:

There exist $\ell \in \bbN$ and $C \ge 0$ such that $\expe[T^\ell] < \infty$ and for all $n \in \bbZ^+$, $|Y_n| \le C \cdot (n+1)^{\ell}$ almost surely.

By $\expe[T^\ell] < \infty$ where $\ell \ge 1$, we know that $\prob[T < \infty] = 1$.

Then similar to the proof of <ref>, we know that $Y_T$ is well-defined almost surely and $\prob[\lim_{n \to \infty} Y_n = Y_T] = 1$.

On the other hand, we have

\begin{align*}
|Y_n| & = |Y_{\min(T,n)}| \le C \cdot (\min(T,n)+1)^\ell \le C \cdot (T + 1)^\ell, \mathrm{a.s.}
\end{align*}

Recall that $\expe[T^\ell] < \infty$.

Then $\expe[(T+1)^\ell] = \expe[T^\ell + O(T^{\ell-1})] < \infty$.

By <ref>, with the function $g$ set to $\lambda\omega. C \cdot (T(\omega)+1)^\ell$, we know that $\lim_{n \to \infty} \expe[Y_n] = \expe[Y_T]$.

By <ref>, we have $\expe[Y_n] = \expe[Y_0]$ for all $n\in\bbZ^+$ thus we conclude that $\expe[Y_0] = \expe[Y_T]$.

§ TRACE-BASED REASONING ON MOMENTS

We start with the fundamental composition property for moment semirings.

For all $u,v \in \calR$, it holds that

\[
\begin{split}
&  \tuple{ 1, (u + v), (u + v)^2, \cdots, (u + v)^m } \\
= {} & \tuple{ 1, u, u^2, \cdots,u^m) \otimes (1, v, v^2, \cdots, v^m },
\end{split}
\]

where $u^n$ is an abbreviation for $\prod_{i=1}^n u$ for $n \in \bbZ^+, u \in \calR$.

Observe that

\begin{align*}
RHS_k & = \sum_{i=0}^k \binom{k}{i} \times (u^i \cdot v^{k-i}).
\end{align*}

We prove by induction on $k$ that $(u + v)^k = RHS_k$.

* $k=0$: Then $(u + v)^0 = 1$.

On the other hand, we have

\[
RHS_0 = \binom{0}{0} \times (u^0 \cdot v^0) = 1 \times (1 \cdot 1) = 1.
\]

* Suppose that $(u + v)^k = RHS_k$.

\begin{align*}
& (u + v)^{k+1} \\
={} & (u + v) \cdot (u + v)^k \\
={} & (u + v) \cdot \sum_{i=0}^k \binom{k}{i} \times (u^i \cdot v^{k-i}) \\
={} & \sum_{i=0}^k \binom{k}{i} \times (u^{i+1} \cdot v^{k-i}) +  \sum_{i=0}^k \binom{k}{i} \times (u^i \cdot v^{k-i+1}) \\
={} & \sum_{i=1}^{k+1} \binom{k}{i-1} \times (u^i \cdot v^{k-i+1}) + \sum_{i=0}^k \binom{k}{i} \times (u^i \cdot v^{k-i+1}) \\
={} & \sum_{i=0}^{k+1} (\binom{k}{i-1} + \binom{k}{i}) \times (u^i \cdot v^{k-i+1}) \\
={} & \sum_{i=0}^{k+1} \binom{k+1}{i} \times (u^k \cdot v^{k-i+1}) \\
={} & RHS_{k+1}.
\end{align*}

We also show that $\otimes$ and $\oplus$ are monotone if the operations of the underlying semiring are monotone.

Let $\calR=(|\calR|,{\le},{+},{\cdot},0,1)$ be a partially ordered semiring.

If $+$ and $\cdot$ are monotone with respect to $\le$, then $\otimes$ and $\oplus$ in the moment semiring $\calM_\calR^{(m)}$ are also monotone with respect to $\aord$.

It is straightforward to show $\oplus$ is monotone.

For the rest of the proof, without loss of generality, we show that $\many{\tuple{u_k}_{0 \le k \le m}} \otimes \many{\tuple{v_k}_{0 \le k \le m}} \aord \many{\tuple{u_k}_{0 \le k \le m}} \otimes \many{\tuple{w_k}_{0 \le k \le m}}$ if $\many{\tuple{v_k}_{0 \le k \le m}} \aord \many{\tuple{w_k}_{0 \le k \le m}}$.

By the definition of $\aord$, we know that $v_k \le w_k$ for all $k=0,1,\cdots,m$.

Then for each $k$, we have

\begin{align*}
& (\many{\tuple{u_k}_{0 \le k \le m}} \otimes \many{\tuple{v_k}_{0 \le k \le m}})_k \\
={} & \sum_{i=0}^k \binom{k}{i} \times (u_i \cdot v_{k - i}) \\
{}\le{} & \sum_{i=0}^k \binom{k}{i} \times (u_i \cdot w_{k - i}) \\
={} & (\many{\tuple{u_k}_{0 \le k \le m}} \otimes \many{\tuple{w_k}_{0 \le k \le m}})_k.
\end{align*}

Then we conclude by the definition of $\aord$.

As we allow potential functions to be interval-valued, we show that the interval semiring $\calI$ satisfies the monotonicity required in <ref>.

The operations $+_\calI$ and $\cdot_\calI$ are monotone with respect to $\le_\calI$.

It is straightforward to show $+_\calI$ is monotone.

For the rest of the proof, it suffices to show that $[a,b] \cdot_\calI [c,d] \le_\calI [a',b'] \cdot_\calI [c,d]$ if $[a,b] \le_\calI [a',b']$, i.e., $[a,b] \subseteq [a',b']$ or $a \ge a', b\le b'$.

We claim that $\min S_{a,b,c,d} \ge \min S_{a',b',c,d}$, i.e., $\min \{ ac, ad, bc, bd\} \ge \min \{ a'c,a'd, b'c,b'd \}$.

* If $0 \le c \le d$: Then $ac \le bc$, $ad \le bd$, $a'c \le b'c$, $a'd \le b' d$.

It then suffices to show that $\min \{ ac, ad\} \ge \min\{ a'c, a'd\}$.

Because $d \ge c \ge 0$ and $a \ge a'$, we conclude that $ac \ge a'c$ and $ad \ge a'd$.

* If $c < 0 \le d$: Then $ac \ge bc$, $ad \le bd$, $a'c \ge b'c$, $a'd \ge b'd$.

It then suffices to show that $\min \{ bc,ad\} \ge \min\{ b'c,a'd\}$.

Because $d \ge 0 > c$ and $a \ge a', b \le b'$, we conclude that $bc \ge b'c$ and $ad \le a'd$.

* If $c \le d < 0$: Then $ac \ge bc$, $ad \ge bd$, $a'c \ge b'c$, $a'd \ge b'd$.

It then suffices to show that $\min\{ bc,bd\} \ge \min\{b'c,b'd\}$.

Because $0 > d \ge c$ and $b \le b'$, we conclude that $bc \ge b'c$ and $bd \ge b'd$.

In a similar way, we can also prove that $\max S_{a,b,c,d} \le \max S_{a',b',c,d}$.

Therefore, we show that $\cdot_\calI$ is monotone.

If $\{[a_n,b_n]\}_{n \in \bbZ^+}$ is a montone sequence in $\calI$, i.e., $[a_0,b_0] \le_\calI [a_1,b_1] \le_\calI \cdots \le_\calI [a_n,b_n] \le_\calI \cdots$, and $[a_n,b_n] \le_\calI [c,d]$ for all $n \in \bbZ^+$.

Let $[a,b] = \lim_{n \to \infty} [a_n,b_n]$ (the limit is well-defined by the monotone convergence theorem for series).

Then $[a,b] \le_\calI [c,d]$.

By the definition of $\le_\calI$, we know that $\{a_n\}_{n \in \bbZ^+}$ is non-increasing and $\{b_n\}_{n \in \bbZ^+}$ is non-decreasing.

Because $a_n \ge c$ for all $n \in \bbZ^+$, we conclude that $\lim_{n \to \infty} a_n \ge c$.

Because $b_n \le d$ for all $n \in \bbZ^+$, we conclude that $\lim_{n \to \infty} b_n \le d$.

Thus we conclude that $[a,b] \le_\calI [c,d]$.

Recall that we extend the notions of $A_n,\Phi_n,Y_n$ with intervals for higher moments as follows:

\begin{align*}
A_n(\omega) & \defeq \many{\tuple{[\alpha_n^k,\alpha_n^k]}_{0 \le k \le m}} ~ \text{where}~\omega_n =\tuple{\_,\_,\_,\alpha_n}, \\
\Phi_n(\omega) & \defeq \phi(\omega_n), \\
Y_n(\omega) & \defeq A_n(\omega) \otimes \Phi_n(\omega)  .
\end{align*}

Note that in the definition of $Y_n$, we use $\otimes$ to compose the powers of the accumulated cost at step $n$ and the potential function that stands for the moments of the accumulated cost of the rest of the computation.

We now extend some of the previous results on first moments to higher moments with intervals.

Let $X,Y : \Omega \to \calM_\calI^{(m)}$ be integrable.

Then $\expe[X \oplus Y] = \expe[X] \oplus \expe[Y]$.

Appeal to linearity of expectations and the fact that $\oplus$ is the pointwise extension of $+_\calI$, as well as $+_\calI$ is the pointwise extension of numeric addition.

If $X:\Omega \to \calM_\calI^{(m)}$ is $\calG$-measurable and bounded, a.s., as well as $X(\omega)=\many{\tuple{[a_k(\omega),a_k(\omega)]}_{0 \le k \le m}}$ for all $\omega \in \Omega$, then $\expe[X \otimes Y \mid \calG] = X \otimes \expe[Y \mid \calG]$ almost surely.

Fix $\omega \in \Omega$. Let $Y(\omega) = \many{\tuple{[b_k(\omega),c_k(\omega)]}_{0 \le k \le m}}$.

Then we have

\begin{align*}
& \expe[X \otimes Y \mid \calG](\omega) \\
={} & \expe[\many{\tuple{[a_k,a_k]}_{0 \le k \le m}} \otimes \many{\tuple{[b_k,c_k]}_{0 \le k \le m}} \mid \calG](\omega) \\
={} & \expe[\many{\tuple{ \textstyle \mathop{\sum_\calI}_{i=0}^k \binom{k}{i} \times_\calI ([a_i,a_i] \cdot_\calI [b_{k-i},c_{k-i}])}_{0 \le k \le m} } \mid \calG](\omega) \\
={} & \many{ \tuple{ \expe[ \textstyle \mathop{\sum_\calI}_{i=0}^k \binom{k}{i} \times_\calI ([a_i,a_i] \cdot_\calI [b_{k-i},c_{k-i}]) \mid \calG](\omega)  }_{0 \le k \le m}} \\
={} & \many{\tuple{ \textstyle \mathop{\sum_\calI}_{i=0}^k \binom{k}{i} \times_\calI  \expe[[a_i,a_i] \cdot_\calI [b_{k-i},c_{k-i}] \mid \calG](\omega) }_{0 \le k \le m} }.
\end{align*}

On the other hand, we have

\begin{align*}
& X(\omega) \otimes \expe[Y \mid \calG](\omega) \\
={} & \tuple{X_0(\omega),\cdots,X_m(\omega)} \otimes \expe[\tuple{Y_0,\cdots,Y_m} \mid \calG](\omega) \\
={} & \tuple{X_0(\omega),\cdots,X_m(\omega)} \otimes \tuple{\expe[Y_0 \mid \calG](\omega), \cdots,\expe[Y_m\mid\calG](\omega)} \\
={} & \many{ \tuple{ \textstyle \mathop{\sum_\calI}_{i=0}^k \binom{k}{i} \times_\calI (X_i(\omega) \cdot_\calI \expe[Y_{k-i} \mid \calG](\omega) ) }_{0 \le k \le m} } \\
={} & \many{ \tuple{ \textstyle \mathop{\sum_\calI}_{i=0}^k \binom{k}{i} \! \times_\calI \! ([a_i(\omega),a_i(\omega)] \! \cdot_\calI \! \expe[[b_{k-i},c_{k-i}] \mid \calG](\omega) ) }_{0 \le k \le m} }.
\end{align*}

Thus, it suffices to show that for each $i$, $\expe[[a_i,a_i] \cdot_\calI [b_{k-i},c_{k-i}] \mid \calG] = [a_i,a_i] \cdot_\calI \expe[[b_{k-i},c_{k-i}] \mid \calG]$ almost surely.

For $\omega$ such that $a_i(\omega) \ge 0$:

\begin{align*}
& \expe[[a_i,a_i] \cdot_\calI [b_{k-i},c_{k-i}] \mid \calG](\omega) \\
={} & \expe[[a_ib_{k-i}, a_ic_{k-i}] \mid \calG](\omega) \\
={} & [\expe[a_ib_{k-i} \mid \calG](\omega), \expe[a_ic_{k-i} \mid \calG](\omega)] \\
={} & [a_i(\omega) \cdot \expe[b_{k-i} \mid \calG](\omega), a_i(\omega) \cdot \expe[c_{k-i} \mid \calG](\omega)], a.s., \\
={} & [a_i(\omega),a_i(\omega)] \cdot_\calI \expe[[b_{k-i},c_{k-i}] \mid \calG](\omega).
\end{align*}

For $\omega$ such that $a_i(\omega) < 0$:

\begin{align*}
& \expe[[a_i,a_i] \cdot_\calI [b_{k-i},c_{k-i}] \mid \calG](\omega) \\
={} & \expe[[a_ic_{k-i}, a_ib_{k-i}] \mid \calG](\omega) \\
={} & [\expe[a_ic_{k-i} \mid \calG](\omega), \expe[a_ib_{k-i} \mid \calG](\omega)] \\
={} & [a_i(\omega) \cdot \expe[c_{k-i} \mid \calG](\omega), a_i(\omega) \cdot \expe[b_{k-i} \mid \calG](\omega)], a.s., \\
={} & [a_i(\omega),a_i(\omega)] \cdot_\calI \expe[[b_{k-i},c_{k-i}] \mid \calG](\omega).
\end{align*}

For all $n \in \bbZ^+$, it holds that

\[
\expe[Y_{n+1} \mid \calF_n] \aord Y_n, \text{a.s.}
\]

Similar to the proof of <ref>, we know that $\{A_n\}_{n \in \bbZ^+}$ and $\{\Phi_n\}_{n \in \bbZ^+}$ are adapted to $\{\calF_n\}_{n \in \bbZ^+}$.

By the property of the operational semantics, we know that $\alpha_n(\omega) \le C \cdot n$ almost surely for some $C \ge 0$.

Then using <ref>, we have

\begin{align*}
& \expe[Y_{n+1} \mid \calF_n](\omega) \\
= {} &\expe[A_{n+1} \otimes \Phi_{n+1} \mid \calF_n](\omega) \\
={} &  \expe[A_n \! \otimes \! \many{\tuple{[(\alpha_{n+1}\!-\!\alpha_n)^k,(\alpha_{n+1}\!-\!\alpha_n)^k]}_{0 \le k \le m}} \!\otimes\! \Phi_{n+1} \mid \calF_n ](\omega) \\
={} & A_n(\omega) \!\otimes\! \expe[\many{\tuple{[(\alpha_{n+1}\!-\!\alpha_n)^k,(\alpha_{n+1}\!-\!\alpha_n)^k]}_{0 \!\le\! k \!\le\! m}} \!\otimes\! \Phi_{n+1} \!\mid\! \calF_n ](\omega), \\
& a.s., \\
={} & A_n(\omega) \!\otimes\! \expe[\many{\tuple{[(\alpha_{n+1}\!-\!\alpha_n)^k,(\alpha_{n+1}\!-\!\alpha_n)^k]}_{0 \!\le\! k \!\le\! m}} \!\otimes\! \phi(\omega_{n+1}) \!\mid\! \calF_n ].
\end{align*}

Recall the property of the expected-potential function $\phi$ in <ref>.

Then by <ref> with <ref>, we have

\begin{align*}
& \expe[Y_{n+1} \mid \calF_n](\omega) \\
{}\aord{} &  A_n(\omega) \otimes \phi(\omega_n), a.s., \\
={} & A_n(\omega) \otimes \Phi_n(\omega) \\
={} & Y_n.
\end{align*}

As a corollary, we have $\expe[Y_n] \aord \expe[Y_0]$ for all $n \in \bbZ^+$.

Now we prove the following extension of OST to deal with interval-valued potential functions.

Let $\norm{\many{\tuple{[a_k,b_k]}_{0 \le k \le m}}}_\infty \defeq \max_{0 \le k \le m} \{ \max\{ |a_k|,|b_k| \} \} $.

If $\expe[\norm{Y_n}_\infty] < \infty$ for all $n \in \bbZ^+$, then $\expe[Y_T]$ exists and $\expe[Y_T] \aord \expe[Y_0]$ in the following situation:

There exist $\ell \in \bbN$ and $C \ge 0$ such that $\expe[T^\ell] < \infty$ and for all $n \in \bbZ^+$, $\norm{Y_n}_\infty \le C \cdot (n+1)^{\ell}$ almost surely.

By $\expe[T^\ell] < \infty$ where $\ell \ge 1$, we know that $\prob[T < \infty] = 1$.

Then similar to the proof of <ref>, we know that $Y_T$ is well-defined almost surely and $\prob[\lim_{n \to \infty} Y_n = Y_T] = 1$.

On the other hand, $Y_n(\omega)$ can be treated as a vector of real numbers.

Let $a_n : \Omega \to \bbR$ be a real-valued component in $Y_n$.

Because $\expe[\norm{Y_n}_\infty] < \infty$ and $\norm{Y_n}_\infty \le C \cdot (n+1)^\ell$ almost surely, we know that $\expe[|a_n|] \le \expe[\norm{Y_n}_\infty] < \infty$ and $|a_n| \le \norm{Y_n}_\infty \le C \cdot (n+1)^\ell$ almost surely.

\[
|a_n| = |a_{\min(T,n)}| \le C \cdot (\min(T,n) + 1)^\ell \le C \cdot (T+1)^\ell, \text{a.s.}
\]

Recall that $\expe[T^\ell] < \infty$.

Then $\expe[(T+1)^\ell] = \expe[T^\ell + O(T^{\ell-1})] < \infty$.

By <ref>, with the function $g$ set to $\lambda\omega. C \cdot (T(\omega) + 1)^\ell$, we know that $\lim_{n \to \infty} \expe[a_n] = \expe[a_T]$.

Because $a_n$ is an arbitrary real-valued component in $Y_n$, we know that $\lim_{n \to \infty} \expe[Y_n] = \expe[Y_T]$.

By <ref>, we know that $\expe[Y_n] \aord \expe[Y_0]$ for all $n \in \bbZ^+$.

By <ref>, we conclude that $\lim_{n \to \infty} \expe[Y_n] \aord \expe[Y_0]$, i.e., $\expe[Y_T] \aord \expe[Y_0]$.

The classic OST from probability theory has been used to reason about probabilistic programs:

Let $\{X_n\}_{n \in \bbZ^+}$ be a super-martingale (resp., sub-martingale).
If $\expe[|X_n|] < \infty$ for all $n \in \bbZ^+$, then $\expe[X_T]$ exists and $\expe[X_T] \le \expe[X_0]$ (resp., $\expe[X_T] \ge \expe[X_0]$) in each of the following situations:

*  $T$ is almost-surely bounded;

*  $\expe[T]\!<\!\infty$ and for some $C \ge 0$, $\expe[|X_{n+1} -X_n| \mid X_n] \le C$ almost surely, for all $n \in \bbZ^+$;

*  for some $C \ge 0$, $|X_n| \le C$ almost surely, for all $n \in \bbZ^+$.

Note that from <ref> to <ref> in <ref>, the constraint on the stopping time $T$ is getting weaker while the constraint on $\{X_n\}_{n \in \bbZ^+}$ is getting stronger.

However, the classic OST is not suitable for higher-moment analysis:
<ref> cannot handle almost-surely terminating programs;
<ref> relies on bounded difference, which is too restrictive, because the change in second moments can be non-constant even if the change in first moments can be bounded by a universal constant;
and <ref> basically requires the accumulated cost can be bounded by a universal constant, which is also not practical.

Recall the random-walk program in <ref>.