Static Factorisation of Probabilistic Programs With User-Labelled Sample Statements and While Loops

TL;DR

This paper introduces a static analysis method for probabilistic programs with user-labelled sample statements and loops, enabling graphical representation and optimizations for inference algorithms.

Contribution

It extends operational semantics and control-flow graph translation to represent unbounded variable dependencies in probabilistic programs with loops and labels.

Findings

Provides a sound static factorisation for programs with loops and labels.

Enables optimizations that reduce variance and improve inference speed.

Empirically outperforms existing inference techniques.

Abstract

It is commonly known that any Bayesian network can be implemented as a probabilistic program, but the reverse direction is not so clear. In this work, we address the open question to what extent a probabilistic program with user-labelled sample statements and while loops - features found in languages like Gen, Turing, and Pyro - can be represented graphically. To this end, we extend existing operational semantics to support these language features. By translating a program to its control-flow graph, we define a sound static analysis that approximates the dependency structure of the random variables in the program. As a result, we obtain a static factorisation of the implicitly defined program density, which is equivalent to the known Bayesian network factorisation for programs without loops and constant labels, but constitutes a novel graphical representation for programs that define an unbounded number of random variables via loops or dynamic labels. We further develop a sound program slicing technique to leverage this structure to statically enable three well-known optimisations for the considered program class: we reduce the variance of gradient estimates in variational inference and we speed up both single-site Metropolis Hastings and sequential Monte Carlo. These optimisations are proven correct and empirically shown to match or outperform existing techniques.

Tables

Table 1. Table 1. Average runtime of computing a single LMH iteration in microseconds (see Section   5.2 ). The column ”Baseline” reports the runtime when computing p ​ ( 𝚝𝚛 ′ ) p(\mathtt{tr}^{\prime}) from scratch. The column ”Static Analysis” reports the runtime when updating the old density by computing log ⁡ p ​ ( 𝚝𝚛 ′ ) − log ⁡ p ​ ( 𝚝𝚛 ) \log p(\mathtt{tr}^{\prime})-\log p(\mathtt{tr}) with the generated sub-programs (see Section   5.1 ). The column ”Fixed-Finite” reports the runtime of updating the density with a hand-written factorisation that is representative of PPLs that explicitly construct a graphical model. The speed-up in relation to computing from scratch is given in the ”speed-up” and ”sp.up” columns. The number of random variables defined in each model can be found in column ”#RVs”. The number of samples produced and the acceptance rate are reported in the N N and A ​ R AR column, respectively.

				Baseline	Static Analysis			Fixed-Finite
Model	#RVs	$N$	$AR$	$\mu s$	$\mu s$	speed-up		$\mu s$	sp.up
Aircraft	$\infty$	$1e5$	0.47	5.861	4.605	1.29	$\nearrow$	-	-
Bayesian Network	37	$1e5$	0.93	25.580	4.680	5.47	$\uparrow \uparrow$	3.344	7.65
Captcha	$\infty$	$1e3$	0.81	774.287	355.382	2.18	$\uparrow \uparrow$	-	-
Dirichlet process	$\infty$	$1e4$	0.53	89.939	42.664	2.12	$\uparrow \uparrow$	-	-
Geometric	$\infty$	$5e5$	0.62	2.419	2.475	0.97	$\to$	-	-
GMM (fixed #clusters)	209	$5e4$	0.60	71.136	7.027	10.01	$\uparrow \uparrow$	4.564	15.58
GMM (variable #clusters)	$\infty$	$5e4$	0.57	74.709	9.549	7.81	$\uparrow \uparrow$	-	-
HMM	101	$1e5$	0.74	33.729	24.016	1.39	$\nearrow$	2.551	13.22
Hurricane	5	$1e6$	0.86	2.121	1.918	1.13	$\to$	-	-
LDA (fixed #topics)	551	$1e4$	0.81	209.046	24.526	8.52	$\uparrow \uparrow$	18.341	11.40
LDA (variable #topics)	$\infty$	$1e4$	0.83	218.614	25.732	8.50	$\uparrow \uparrow$	-	-
Linear regression	102	$1e5$	0.53	4.477	8.006	0.56	$\downarrow \downarrow$	7.504	0.60
Marsaglia	$\infty$	$5e5$	0.93	1.699	2.159	0.79	$\searrow$	-	-
PCFG	$\infty$	$1e5$	0.75	10.220	8.437	1.21	$\to$	-	-
Pedestrian	$\infty$	$1e5$	0.37	3.709	3.645	1.02	$\to$	-	-
Urn	$\infty$	$1e5$	0.52	11.979	6.354	1.88	$\nearrow$	-	-

Table 2. Table 2. Relative speed-up for the LMH algorithm achieved with our static approach, combinators in Gen, and the C3 algorithm in WebPPL.

Model	ours	Gen	C3	Model	ours	Gen	C3
Aircraft	1.29	1.41	1.02	Hurricane	1.13	0.51	0.93
Bayesian Network	5.47	4.14	0.93	LDA (fixed #topics)	8.52	6.50	4.97
Captcha	2.18	2.93	1.05	LDA (variable #topics)	8.50	5.92	7.26
Dirichlet process	2.12	-	1.24	Linear regression	0.56	0.43	0.44
Geometric	0.97	-	0.74	Marsaglia	0.79	-	0.53
GMM (fixed #clusters)	10.01	6.26	1.90	PCFG	1.21	-	0.79
GMM (variable #clusters)	7.81	5.45	1.60	Pedestrian	1.02	-	0.52
HMM	1.39	1.31	3.82	Urn	1.88	2.05	0.69

Table 3. Table 3. Average variance of the standard and improved ELBO gradient estimator in BBVI. The variance is averaged over all parameters, 1000 iterations, and 100 samples per iteration. The ”cost” column reports the factor by which the runtime increased for the improved estimator. For reference, we include the variance reduction factor achieved by Pyro’s TraceGraph_ELBO .

	Baseline	Static Analysis				Pyro
Model	avg. var.	avg. var.	red. factor		cost	red. factor	cost
Aircraft	7.922e+03	7.969e+01	99.29	$\downarrow \downarrow$	1.89	-	-
Bayesian Network	4.160e+02	4.371e+00	95.06	$\downarrow \downarrow$	1.94	95.09	7.52
Geometric	1.634e-01	1.726e-01	0.94	$\to$	0.53	-	-
GMM (fixed #clusters)	8.043e+09	1.014e+06	8092.80	$\downarrow \downarrow$	2.26	9848.79	8.69
HMM	2.116e+03	1.279e+03	1.65	$\to$	7.49	1254.24	8.67
HMM (no loop)	2.116e+03	2.349e+00	901.01	$\downarrow \downarrow$	1.50	1254.24	8.64
Hurricane	3.362e-01	3.609e-01	0.93	$\to$	0.81	-	-
LDA (fixed #topics)	7.229e+04	1.172e+00	61646.02	$\downarrow \downarrow$	5.92	59986.61	11.00
Linear regression	3.350e+05	3.341e+05	1.00	$\to$	1.14	1.00	14.39
Marsaglia	7.274e-02	1.283e-02	5.38	$\searrow$	0.51	-	-
PCFG	8.838e-03	1.989e-02	0.56	$\to$	0.94	-	-
Pedestrian	3.108e+06	3.111e+06	1.00	$\to$	2.01	-	-

Table 4. Table 4. Average runtime of naive versus iterative implementations of SMC with 100 particles. For reference, we report the speed-up achieved by WebPPL’s CPS transformation.

		Baseline	Static Analysis			WebPPL
Model	#Data	$ms$	$ms$	speed-up		speed-up
Dirichlet process	50	238.323	22.374	10.57	$\uparrow \uparrow$	8.96
GMM (fixed #clusters)	100	166.480	31.260	5.26	$\uparrow \uparrow$	9.23
GMM (variable #clusters)	100	189.529	31.263	5.95	$\uparrow \uparrow$	9.18
HMM	50	31.951	7.394	4.13	$\uparrow \uparrow$	7.66
HMM (no loop)	50	49.649	7.420	6.63	$\uparrow \uparrow$	6.79
LDA (fixed #topics)	262	1319.553	397.074	3.33	$\uparrow \uparrow$	7.25
LDA (variable #topics)	262	1412.182	410.781	3.44	$\uparrow \uparrow$	8.18

Equations

$$\sigma ::=x_1\mapsto V_1,\dots ,x_n\mapsto V_n,x_i\text{distinct},V_i\in \mathcal{V}.$$

$$\sigma (x_i)=V_i,\sigma (c)=c,\sigma (g(E_1,\dots ,E_n))=g(\sigma (E_1),\dots ,\sigma (E_n)).$$

$$\frac{}{(\sigma ,\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierskip}){\Downarrow }^{\mathtt{tr}}\sigma \frac{}{(\sigma ,x=E){\Downarrow }^{\mathtt{tr}}\sigma [x\mapsto \sigma (E)]}\frac{(\sigma ,S_1){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }({\sigma }^{\prime },S_2){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime \prime }}{(\sigma ,S_1;S_2){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime \prime }}$$

$$\frac{\sigma (E)=\mathtt{true}(\sigma ,S_1){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }}{(\sigma ,\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierif}E\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierthen{S}_1\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierelse{S}_2){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }\frac{\sigma (E)=\mathtt{false}(\sigma ,S_2){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }}{(\sigma ,\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierif}E\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierthen{S}_1\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierelse{S}_2){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }$$

$$\frac{\sigma (E)=\mathtt{false}}{(\sigma ,\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierwhile}E\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierdoS){\Downarrow }^{\mathtt{tr}}\sigma \frac{\sigma (E)=\mathtt{true}(\sigma ,(S;\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierwhile}{}E\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierdoS)){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }(\sigma ,\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierwhileE\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierdoS){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }$$

$$(\sigma ,\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierwhileE\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifierdoS){\Downarrow }^{\mathtt{tr}}{\sigma }^{\prime }⟺$$

$$\frac{\forall i:\sigma (E_i)=V_i\wedge V_i\ne \mathtt{null}{V}_0\in \mathbf{Strings}V=\mathtt{tr}(V_0)\wedge V\ne \mathtt{null}}{(\sigma ,x=\text{mbox}\text{leavevmode}\text{lstinline}\text{lst@@@set@language}\text{lst@@@set@numbers}\text{lst@@@set@frame}\text{lst@@@set@rulecolor}\text{@listingGroup}lt{x}_{l}s{t}_{i}dentifiersample}(E_0,f(E_1,\dots ,E_n))){\Downarrow }^{\mathtt{tr}}\sigma [x\mapsto V,\text{p}\mapsto \sigma (\text{p})\times {\text{pdf}}_f(V;V_1,\dots ,V_n)]$$

$${\bar{p}}_S(\mathtt{tr},{\sigma }_0):=\{\begin{array}{ll}\sigma (\text{p})& \text{if }\exists \sigma \in \Sigma :({\sigma }_0,S){\Downarrow }^{\mathtt{tr}}\sigma \\ \text{undefined}& \text{otherwise}\end{array}$$

$$p(\mathtt{tr})$$

$$\times {\text{pdf}}_{\text{Bernoulli}}(\mathtt{tr}(\text{"x"});\mathtt{tr}(\text{"p"}))$$

$$\times \left({\delta }_{\mathtt{tr}(\text{"x"})}{\text{pdf}}_{\text{Bernoulli}}(\mathtt{tr}(\text{"y"});0.25)+(1-{\delta }_{\mathtt{tr}(\text{"x"})})\right)$$

$$\times \left({\delta }_{\mathtt{tr}(\text{"x"})}+(1-{\delta }_{\mathtt{tr}(\text{"x"})}){\text{pdf}}_{\text{Bernoulli}}(\mathtt{tr}(\text{"z"});0.75)\right).$$

$$(\sigma ,N)\stackrel{\mathtt{tr}}{⟶}({\sigma }^{\prime },N^{\prime })$$

$$\frac{N=\text{START}{N}^{\prime }=\text{successor}(N)}{(\sigma ,N)\stackrel{\mathtt{tr}}{⟶}(\sigma ,N^{\prime })}\frac{N=\text{Assign}(x=E)N^{\prime }=\text{successor}(N)}{(\sigma ,N)\stackrel{\mathtt{tr}}{⟶}(\sigma [x\mapsto \sigma (E)],N^{\prime })}$$

$$\frac{N=\text{Join}{N}^{\prime }=\text{successor}(N)}{(\sigma ,N)\stackrel{\mathtt{tr}}{⟶}(\sigma ,N^{\prime })}\frac{N=\text{Branch}(E)\sigma (E)=\mathtt{true}}{(\sigma ,N)\stackrel{\mathtt{tr}}{⟶}(\sigma ,\text{consequent}(N))}\frac{N=\text{Branch}(E)\sigma (E)=\mathtt{false}}{(\sigma ,N)\stackrel{\mathtt{tr}}{⟶}(\sigma ,\text{alternative}(N))}$$

$${\frac{\begin{aligned} {N=\text{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},f(E_{1},\dots,E_{n})))\quad N^{\prime}=\text{successor}(N)}\\ {\forall i\colon\sigma(E_{i})=V_{i}\land V_{i}\neq\mathtt{null}\quad V_{0}\in\mathbf{Strings}\quad V=\mathtt{tr}(V_{0})\land V\neq\mathtt{null}}\end{aligned}}{(\sigma,N)\overset{\mathtt{tr}}{\longrightarrow}(\sigma[x\mapsto V,\textbf{p}\mapsto\sigma(\textbf{p})\times\textnormal{pdf}_{f}(V;V_{1},\dots,V_{n})],N^{\prime})}$$

$${\bar{p}}_G(\mathtt{tr},{\sigma }_0):=\{\begin{array}{ll}\sigma (\text{p})& \text{if}\exists \sigma \in \Sigma ,n\in \mathbb{N}:\forall i=1,\dots ,n:\exists {\sigma }_i,N_i:\\ & ({\sigma }_0,\text{START})\stackrel{\mathtt{tr}}{⟶}\dots \stackrel{\mathtt{tr}}{⟶}({\sigma }_i,N_i)\stackrel{\mathtt{tr}}{⟶}\dots \stackrel{\mathtt{tr}}{⟶}(\sigma ,\text{END})\\ \text{undefined}& \text{otherwise.}\end{array}$$

$${\bar{p}}_S(\mathtt{tr},{\sigma }_0)={\bar{p}}_G(\mathtt{tr},{\sigma }_0).$$

$$\text{RD}(N,x)=\{N^{\prime }\in \text{AssignNodes}(G,x):\exists \text{ path }(N^{\prime },N_1,\dots ,N_n,N),N_i\notin \text{AssignNodes}(G,x)\}.$$

$$\text{BP}(N)=\{B:(B,J)\in \text{BranchJoin}(G)\wedge \exists \text{ path }(B,\dots ,N,\dots ,J)\}.$$

$$\text{addresses}(N)=\{\sigma (E_0):\sigma \in \Sigma \}\subseteq \mathbf{Strings}.$$

$$\text{prov}(N,E)={\bigcup }_{y\in \text{vars}(E)}\text{prov}(N,y).$$

$$N_7=\text{Assign}(\mathtt{x}=\dots ),V_{\mathtt{s}}^{N_7}(\mathtt{tr})=\mathtt{tr}(\text{"s"}),$$

$$V_{\mathtt{m}}^{N_7}(\mathtt{tr})=\{\begin{array}{ll}\mathtt{tr}(\text{"mu"})& \text{if }\mathtt{tr}(\text{"b"})=\mathtt{true},\\ 1& \text{otherwise.}\end{array}$$

$$\forall \alpha \in A:\mathtt{tr}(\alpha )={\mathtt{tr}}^{\prime }(\alpha )⟹f(\mathtt{tr})=f({\mathtt{tr}}^{\prime }).$$

$$A\subseteq A^{\prime }⟹[\mathcal{T}{|}_A\to B]\subseteq [\mathcal{T}{|}_{A^{\prime }}\to B].$$

$$\mathtt{tr}\mapsto h(f_1(\mathtt{tr}),f_2(\mathtt{tr}))\in [\mathcal{T}{|}_{A_1\cup A_2}\to C].$$

$$V_x^N\in [\mathcal{T}{|}_{\text{prov}(N,x)}\to \mathcal{V}],$$

$$({\sigma }_0,\text{START})\stackrel{\mathtt{tr}}{⟶}\cdots \stackrel{\mathtt{tr}}{⟶}({\sigma }_i,N_i)\stackrel{\mathtt{tr}}{⟶}\cdots \stackrel{\mathtt{tr}}{⟶}({\sigma }_l,\text{END})$$

$${\sigma }_i(x)=V_x^{N_i}(\mathtt{tr}).$$

$$V_c^N(\mathtt{tr})=c,V_{g(E_1,\dots ,E_n)}^N(\mathtt{tr})=g(V_{E_1}^N(\mathtt{tr}),\dots ,V_{E_n}^N(\mathtt{tr})).$$

$$A_k=\text{addresses}(N_k)\cup {\bigcup }_{i=0}^{n_k}\text{prov}(N_k,E_i^k)\cup {\bigcup }_{N^{\prime }\in \text{BP}(N_k)}\text{prov}(N^{\prime },\text{condexp}(N^{\prime })).$$

Full text

Probabilistic programming, operational semantics, static program analysis, factorisation, Bayesian networks

††copyright: acmlicensed††journalyear: 2025††doi: XXXXXXX.XXXXXXX††conference: ACM Conference; July 2025; Washington, DC, USA††isbn: 978-1-4503-XXXX-X/18/06††ccs: Theory of computation Operational semantics††ccs: Theory of computation Program analysis††ccs: Mathematics of computing Bayesian computation††ccs: Mathematics of computing Bayesian networks

1. Introduction

Probabilistic programming provides an intuitive means to specify probabilistic models as programs. Typically, probabilistic programming systems also come with machinery for automatic posterior inference. This enables practitioners without substantial knowledge in Bayesian inference techniques to work with complex probabilistic models.

Besides these practical aspects, there is also a long history of studying probabilistic programming languages (PPLs) from a mathematical perspective (Kozen, 1979; Barthe et al., 2020). It is commonly known that any Bayesian network (Koller and Friedman, 2009) can be implemented as a probabilistic program (Gordon et al., 2014), but the reverse direction is not so clear. The understanding is that first-order probabilistic programs without loops can be compiled to Bayesian networks, because the number of executed sample statements is fixed and finite, such that each sample statement corresponds to exactly one random variable (van de Meent et al., 2018; Borgström et al., 2011; Sampson et al., 2014). Generally, a graphical representation of a probabilistic program is desirable as it lends itself for optimising inference algorithms by exploiting the dependency structure between random variables. In fact, many probabilistic programming languages require the user to explicitly construct the model in a graphical structure (Salvatier et al., 2016; Minka, 2012; McCallum et al., 2009; Lunn et al., 2000; Milch et al., 2007).

However, in the universal probabilistic programming paradigm probabilistic constructs are embedded in a Turing-complete language (Goodman et al., 2012). While loops and recursive function calls easily lead to an unbounded number of executed sample statements which makes the typical compilation to a graph with one node per random variable impossible. Even more, many PPLs allow the user to label each sample statement with a dynamic address. This further complicates matters, because in this case even a single sample statement may correspond to a potentially infinite number of random variables. As a consequence, many approaches resort to dynamic methods to leverage the dependency structure of this kind of probabilistic programs for efficient inference (Wu et al., 2016; Yang et al., 2014; Chen et al., 2014).

Sample statements with dynamic address expressions constitute a valuable language feature, as they allow the definition of random variables whose identifiers differ from the program variables to which values are assigned to. This permits the specification of a broader model class relative to PPLs that lack this feature. Below we give examples of such user-labelled sample statements in Gen (Cusumano-Towner et al., 2019), Pyro (Bingham et al., 2019), and Turing (Ge et al., 2018), where the address expressions are highlighted with red:

⬇

x = {:x => i} ~ Normal(0.,1.) # Gen

x = pyro.sample(f"x{i}", Normal(0.,1.)) # Pyro

x[i] ~ Normal(0.,1.) # Turing

x = sample("x"+str(i), Normal(0.,1.)) # our formal syntax

In this work, we address the open question to what extent a probabilistic program with user-labelled sample statements and while loops can be represented graphically. We prove that that the density of a probabilistic program with $K$ user-labelled sample statements decomposes into $K$ factors. This finite factorisation holds even if the program contains while loops. Moreover, the dependency structure implied by the factorisation can be computed statically which, as we will show, enables new opportunities for the optimisation of inference algorithms.

1.1. Overview

Section 2**. Address-Based Semantics For Probabilistic Programs** generalises the operational semantics of Core Stan (Gorinova et al., 2019). Instead of interpreting the program density $p$ as a function over only a fixed set of finite variables, we extend this approach to interpret the density as a function over traces. In contrast to Gorinova et al. (2019), this allows us to consider programs with dynamic addresses and while loops. These language features greatly expand the class of expressible probabilistic models, for instance, models with an infinite number of random variables. The described semantics can be seen as an operational analog to the denotational models of Pyro (Lee et al., 2019) and Gen (Lew et al., 2023).

Section 3**. Static Factorisation for Programs Without Loops** examines the factorisation of the program density for the simplified case where programs do not contain while loops.

In Section 3.2, building on standard data and control flow analysis techniques, we present a static provenance analysis that finds the dependencies of a program variable $x$ at any program state $\sigma$. The dependencies are described by specifying for which addresses a trace $\mathtt{tr}$ needs to be known to also know the value of $x$ in $\sigma$. This analysis forms the basis for finding the factorisation over addresses and operates on the control-flow graph (CFG) of a program.

In Section 3.3, we derive a novel factorisation from the static analysis, $p(\mathtt{tr})=\prod_{k=1}^{K}p_{k}(\mathtt{tr})$, for programs without loops by explicitly constructing one factor $p_{k}$ per sample statement. If we further restrict the program to constant addresses, we arrive at the language of Gorinova et al. (2019), where our factorisation theorem implies a known equivalence to Bayesian networks. However, we show that already in the case of dynamic addresses, programs are only equivalent to Markov networks, because of the possibility of cyclic dependencies and an infinite number of random variables.

Section 4**. Static Factorisation for Programs With Loops** describes a way to generalise the factorisation result of Section 3 to the much more challenging case where programs may contain while loops. The idea is to mathematically unroll while loops to obtain a directed acyclic graph, the unrolled CFG, that essentially represents all possible program paths. Our main theoretical contribution is to show that even in this case, the density still factorises into a finite number of terms – one factor for each of the $K$ sample statements. Like dynamic addresses, while loops may also be used to model an infinite number of random variables and thus, the resulting factorisation in general implies the equivalence of programs to Markov networks.

Section 5**. Optimising Posterior Inference Algorithms** presents a novel and sound program slicing method building on the factorisation theory established in Section 4. By generating sub-programs for each factor, we can 1. decrease the runtime of a Metropolis Hastings algorithm by speeding up the acceptance probability computation; 2. construct a lower-variance gradient estimator for Black-Box Variational Inference; 3. facilitate an iterative implementation of sequential Monte Carlo. Each of the considered optimisations are known, but have been either realised dynamically or were reserved for PPLs that require the explicit graphical construction of probabilistic models. Our approach is the first to enable them statically for the considered program class. We evaluate our approach on a benchmark set of 16 probabilistic programs, the majority of which make use of while loops to define an unbounded number of random variables.

1.2. Contributions

We briefly summarise the key contributions of this work:

• •

A static provenance analysis in the context of our formal PPL building on standard data and control flow techniques,

• •

A novel density factorisation for the considered program class, which can be computed statically and subsumes a known equivalence to Bayesian networks in a restricted setting,

• •

A novel and sound program slicing method derived from the factorisation theory,

• •

An application of the slicing method to enable well-known optimisations statically and provably correct for the considered program class and three inference algorithms: Metropolis Hastings, variational inference, and sequential Monte Carlo,

• •

An open-source implementation and empirical evaluation of our proposed approach (Böck, 2025).

2. Address-Based Semantics For Probabilistic Programs

First, we introduce the syntax of our formal probabilistic programming language in terms of expressions and statements similar to Gorinova et al. (2019) and Hur et al. (2015). We assume that the constants of this language range over a set of values including booleans, integers, real numbers, immutable real-valued vectors, finite-length strings, and a special $\mathtt{null}$ value. We define $\mathcal{V}$ to be the set of all values. A program contains a finite set of program variables denoted by $x$, $x_{i}$, $y$, or $z$. We further assume a set of built-in functions $g$. We assume that functions $g$ are total and return $\mathtt{null}$ for erroneous inputs. The only non-standard construct in the language are sample statements. For sample statements, the expression $E_{0}$ corresponds to the user-defined sample address, which is assumed to evaluate to a string. The symbol $f$ ranges over a set of built-in distributions which are parameterised by arguments $E_{1},\dots,E_{n}$.

[TABLE]

[TABLE]

2.1. Operational Semantics

The meaning of a program is the density implicitly defined by its sample statements. This density is evaluated for a program trace and returns a real number. A program trace is a mapping from finite-length addresses to numeric values, $\mathtt{tr}\colon\mathbf{Strings}\to\mathcal{V}$. We denote the set of all traces with $\mathcal{T}$. How the value at an address is determined is not important for this work and in practice typically depends on the algorithm that is used to perform posterior inference for the program. For instance, if an address $\alpha$ corresponds to an observed variable, the value $\mathtt{tr}(\alpha)$ is fixed to the data. If an address $\alpha$ does not appear in a program, then its value can be assumed to be $\mathtt{tr}(\alpha)=\mathtt{null}$.

For each trace $\mathtt{tr}$, the operational semantics of our language are defined by the big-step relation $(\sigma,S)\Downarrow^{\mathtt{tr}}\sigma^{\prime}$, where $S$ is a statement and $\sigma\in\Sigma$ is a program state – a finite map from program variables to values

[TABLE]

A state can be naturally lifted to a map from expressions to values $\sigma\colon\mathbf{Exprs}\to\mathcal{V}$:

[TABLE]

The operational semantics for the non-probabilistic part of our language are standard:

[TABLE]

[TABLE]

[TABLE]

Above, $\sigma[x\mapsto V]$ denotes the updated program state $\sigma^{\prime}$, where $\sigma^{\prime}(x)=V$ and $\sigma^{\prime}(y)=\sigma(y)$ for all other variables $y\neq x$. The inference rules are to be interpreted inductively such that the semantics of a while loop can be rewritten as

[TABLE]

where $\textnormal{repeat}_{n}(S)=(S;...;S)$ repeated $n$-times. In particular, if a while loop does not terminate for state $\sigma$ then ${{\nexists\sigma^{\prime}\colon(\sigma,\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{while}}}}}}\;E\;\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{do}}}}}}\;S)\Downarrow^{\mathtt{tr}}\sigma^{\prime}$. This makes the presented operational semantics partial, where errors and non-termination are modelled implicitly. Formally, $\Downarrow^{\mathtt{tr}}$ is the least-relation closed under the stated rules.

At sample statements, we inject the value of the trace at address $V_{0}={\sigma(E_{0})}$. We multiply the density, represented by the reserved variable p, with the value of function $\textnormal{pdf}_{f}$ evaluated at $\mathtt{tr}(V_{0})$. Note that the function $\textnormal{pdf}_{f}$ may be an arbitrary function to real values. However, in the context of probabilistic programming, we interpret $\textnormal{pdf}_{f}$ as the density function of distribution $f$ with values in $\mathbb{R}_{\geq 0}$. Lastly, the value $\mathtt{tr}(V_{0})$ is stored in the program state at variable $x$. The semantics of sample statements are summarised in the following rule:

[TABLE]

Finally, the meaning of a program $S$ is defined as the mapping of trace to density. This definition is well-defined, because the semantics are deterministic for each trace. If a trace leads to errors or non-termination, then the density is undefined.

Definition 0.

The density of a program $S$ is the function $\bar{p}_{S}\colon\mathcal{T}\times\Sigma\to\mathbb{R}_{\geq 0}$ given by

[TABLE]

For the initial program state $\sigma_{0}=(x_{1}\mapsto\mathtt{null},\dots,x_{n}\mapsto\mathtt{null},\textbf{p}\mapsto 1)$, mapping all program variables to $\mathtt{null}$, we define $p_{S}\colon\mathcal{T}\to\mathbb{R}_{\geq 0}$ with $p_{S}(\mathtt{tr})=\bar{p}_{S}(\mathtt{tr},\sigma_{0}).$

Consider the simple probabilistic program in LABEL:lst:semantics (Note that we slightly deviate from the formal syntax in listings for better readability). The density defined by our operational semantics is given below, where $\delta_{v}$ denotes the delta function whose value is 1 if $v$ equals $\mathtt{true}$ and 0 otherwise.

[TABLE]

As the addresses in LABEL:lst:semantics are constant and identical to the program variables, the program is well-described by many existing PPL semantics. However, in LABEL:lst:semantics2, we rewrite the program with dynamic addresses which results in the same above density. Most formal PPLs studied in prior work do not support user-defined addresses for sample statements and cannot interpret a program with three sample statements as a four-dimensional joint density. In contrast, address-based semantics like those of Lee et al. (2019) and Lew et al. (2023), denotationally modelling Pyro and Gen, respectively, correctly interpret the program with dynamic addresses. This is also true for our extensions of Gorinova et al. (2019). While closely related to the aforementioned denotational models, we opted for an operational approach to more directly model the implementations of imperative PPLs.

Note that probabilistic programs with dynamic addresses and while loops pose a significant theoretical challenge, because both cases easily lead to an unbounded number of addresses as can be seen in LABEL:lst:poisson and LABEL:lst:geometric.

2.2. Remark on a Measure-theoretic Interpretation

Even though it suffices to consider the explicit computation of the density of a probabilistic program in this work, a measure-theoretic interpretation of the presented operational semantics is still worth exploring. The density of a measure $\mu$ with respect to some reference measure $\nu$ is formally given by the Radon-Nikodym derivative $p$ such that $\mu(A)=\int_{A}p\;d\nu$. Thus, the measure-theoretic interpretation of the operational semantics is described by a construction of a measure space of traces $(\mathcal{T},\Sigma_{\mathcal{T}},\nu)$ such that the density $p_{S}$ is measurable. Then, $\mu_{S}(A)=\int_{A}p_{S}\;d\nu$ is the measure denoted by the program $S$. As our results do not rely on measure theoretic assumptions, we refer to Appendix C for a sketch of this measure space.

2.3. Remark on Conditioning

In the presented semantics conditioning is only supported via soft-constraints (Staton et al., 2016). This means that the input trace $\mathtt{tr}$ is assumed to be fixed to observed data at certain addresses. If those addresses are encountered during execution at sample statements the data likelihood is multiplied to $\sigma(\textbf{p})$. Informally, the program density $p_{S}$ corresponds to the unnormalised posterior density. Again, a rigorous treatment of the normalisation of $p_{S}$ is not required for our results.

3. Static Factorisation for Programs Without Loops

The goal of this work is to statically find a factorisation of the density implicitly defined by a probabilistic program in terms of addresses. This static analysis will operate on the control-flow graph (CFG) of a program. We begin by considering programs without loops first and describe how we can formally translate such a program to its corresponding CFG. We will equip this CFG with semantics equivalent to the operational semantics of the original program. The following construction is fairly standard but necessary for formally verifying the correctness of the factorisation.

3.1. Control-Flow Graph

A CFG has five types of nodes: start, end, assign, branch, and join nodes. Branch nodes have two successor nodes, end nodes have no successor, all other nodes have one. Only join nodes have multiple predecessors. The CFG contains exactly one start node from which every other node is reachable. Further, it contains exactly one end node that is reachable from any other node. We describe the recursive translation rules for program $S$ as diagrams, which serve to convey the general idea of the construction. For a formal description see Appendix A. Sub-graphs are drawn by circular nodes. Branch and join nodes are drawn as diamond nodes, the latter filled black. Assign, start, and end nodes have rectangular shape. In text, branch nodes are written as $\text{Branch}(E)$. Assign nodes are denoted with $\text{Assign}(x=E)$ and ${\text{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},\dots))$.

• •

The CFG of a skip-statement, ${S=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{skip}}}}}}$, consists of only start and end node:

• •

The CFG of an assignment, $x=E$, or sample statement ${x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},f(E_{1},\dots,E_{n}))$, is a sequence of start, assign, and end node:

• •

The CFG for a sequence of statements, $S=S_{1};S_{2}$, is recursively defined by ”stitching together” the CFGs $G_{1}$ and $G_{2}$ of $S_{1}$ and $S_{2}$:

• •

The CFG of an if statement, ${{{S=(\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{if}}}}}}\;E\;\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{then}}}}}}\;S_{1}\,\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{else}}}}}}\;S_{2})$, is also defined in terms of the CFGs of $S_{1}$ and $S_{2}$ together with a branch node $B=\textnormal{Branch}(E)$ and a join node $J$:

We denote the successor nodes of the branch node $B$ as $\textnormal{consequent}(B)$ and $\textnormal{alternative}(B)$. Furthermore, we say the tuple $(B,J)$ is a branch-join pair and define $\text{BranchJoin}(G)$ to be the set of all branch-join pairs of $G$. Lastly, let $\textnormal{condexp}(B)$ denote the conditional expression $E$.

3.1.1. CFG Semantics

For a control-flow graph $G$ we define the small-step semantics for a trace $\mathtt{tr}$ as a relation that models the transition from node to node depending on the current program state:

[TABLE]

The semantics for the CFG of program $S$ are set-up precisely such that they are equivalent to the operational semantics of program $S$. The transition rules depending on node type and program state are given below, where we denote the unique child of non-branch nodes $N$ in $G$ with successor(N):

[TABLE]

[TABLE]

[TABLE]

Finally, the meaning of a CFG is a mapping from trace to density. It is defined for traces $\mathtt{tr}$ if there exists a path of node transitions from start to end node.

Definition 0.

The density of a CFG $G$ is the function $\bar{p}_{G}\colon\mathcal{T}\times\Sigma\to\mathbb{R}_{\geq 0}$ given by

[TABLE]

For initial program state $\sigma_{0}$ as in Definition 1, we define $p_{G}\colon\mathcal{T}\to\mathbb{R}_{\geq 0}$ with $p_{G}(\mathtt{tr})=\bar{p}_{G}(\mathtt{tr},\sigma_{0}).$

Proposition 0.

For all programs $S$ without while loops and corresponding control-flow graph $G$, it holds that for all traces $\mathtt{tr}$ and all initial states $\sigma_{0}$

[TABLE]

Proof.

By structural induction. See Appendix B.1. ∎

3.2. Static Provenance Analysis

The provenance of variable $x$ in CFG node $N$ are all the addresses $\subseteq\mathbf{Strings}$ that contribute to the computation of its value in node $N$. In this section, we define an algorithm that operates on the CFG and statically determines the provenance for any variable $x$ at any node $N$. As this is a static algorithm, it is an over-approximation of the true provenance. In Section 3.2.1, we will prove soundness of the presented algorithm. We show that indeed the value of variable $x$ at node $N$ can be computed from the values in the trace $\mathtt{tr}$ at the addresses which make up the provenance. But first, we introduce some useful concepts which are similar to definitions found in data-flow analysis (Cooper and Torczon, 2011; Hennessy and Patterson, 2011).

Definition 0.

For a CFG $G$, let $\textnormal{AssignNodes}(G,x)$ be all nodes $N^{\prime}$ in $G$ that assign variable $x$, $N^{\prime}=\textnormal{Assign}(x=\dots)$. The set of reaching definitions for variable $x$ in node $N$ is defined as

[TABLE]

Intuitively, the set of reaching definitions are all nodes that could have written the value of $x$ in the program state before executing node $N$.

Definition 0.

For CFG $G$ and node $N$ the set of branch parents is defined as

[TABLE]

A branch node $B$ and join node $J$ are a branch-join pair $(B,J)\in\textnormal{BranchJoin}(G)$ if they belong to the same if-statement (see the translation rule for if-statements to CFGs).

The set of branch parents are all parent branch nodes of $N$ that determine if the program branch that $N$ belongs to is executed.

Lastly, we define the set of all strings that the address expression of a sample statement may evaluate to.

Definition 0.

For a sample CFG node ${N=\textnormal{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},f(E_{1},\dots,E_{n})))$, we define the set of all possible addresses as

[TABLE]

Now we are able to define the algorithm to statically approximate the provenance of any variable $x$ at any node $N$ denoted by $\textnormal{prov}(N,x)$ in Algorithm 1. This algorithm is inspired by standard methods for dependence analysis. In the algorithm, we use $\textnormal{vars}(E)$, which is defined as the set of all program variables in expression $E$. It is often useful to determine the provenance of an entire expression instead of only a single variable. For expressions $E$, we lift the definition of prov

[TABLE]

The algorithm works as follows. If we want to find $\textnormal{prov}(N,x)$, we first have to find all reaching definitions for $x$. If the reaching definition $N^{\prime}$ is a sample node, ${x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},\dots)$, then the value of $x$ can only depend on the addresses generated by $E_{0}$, $\textnormal{addresses}(N^{\prime})$, and the provenance of $E_{0}$. If $N^{\prime}$ is an assignment node, $x=E$, we have to recursively find the provenance of the expression $E$. Lastly, there can be multiple reaching definitions for $x$ in different branches and the value of $x$ also depends on the branching condition. Thus, we also find all branch parents and recursively determine their provenance and add it to the final result.

In practice the potentially infinite set $\textnormal{addresses}(N)$ can be represented by the CFG node $N$ itself. As remarked on line 8, an implementation of Algorithm 1 would return a set of CFG sample nodes rather than a potentially infinite set of strings. If we call this version of the algorithm provnode, we have $\textnormal{prov}(N,x)=\bigcup_{N^{\prime}\in\textnormal{provnode}(N,x)}\textnormal{addresses}(N^{\prime})$.

We get a clearer sense of the introduced definitions and Algorithm 1 by considering a simple program in Figure 1. The right panel shows the CFG with reaching definitions and branch parents for several nodes drawn as arrows. These arrows illustrate how Algorithm 1 computes the provenance of variables m and s at the sample node for address ”x”. Note that if we were to change the if statement to (if (b == 1) then (m = 1) else (m = 1)), then the address ”b” would still be in the provenance set of m. This illustrates that the presented algorithm can only over-approximate the provenance and not compute it exactly.

3.2.1. Soundness

Recall that the purpose of finding the provenance of variable $x$ at node $N$ is to determine which addresses contribute to the computation of its value. To prove that Algorithm 1 indeed finds the correct set of addresses, first, we have to precisely define what we mean that a set of addresses contribute to the computation. Suppose that there is a function $f$ that maps a trace $\mathtt{tr}$ to a value $v$. The set of addresses $A$ contribute to the computation of value $v$, if any other trace $\mathtt{tr}^{\prime}$ that has the same values as $\mathtt{tr}$ for all addresses $\alpha\in A$, gets mapped to the same value $f(\mathtt{tr}^{\prime})=v$. In other words, if we change the trace $\mathtt{tr}$ at an address $\beta\in\mathbf{Strings}\setminus A$, the value of $f(\mathtt{tr})$ remains unchanged.

Definition 0.

For a function from traces to an arbitrary set $B$, $f:\mathcal{T}\to B$, we write $f\in[\mathcal{T}|_{A}\to B]$ for a subset of addresses $A\subseteq\mathbf{Strings}$, if and only if

[TABLE]

That is, changing the values of $\mathtt{tr}$ at addresses $\beta\in\mathbf{Strings}\setminus A$ does not change the value $f(\mathtt{tr})$.

The next lemma states that we can always over-approximate the set of addresses that contribute to the computation of the values $f(\mathtt{tr})$.

Lemma 0.

If $A\subseteq A^{\prime}$ and $f\in[\mathcal{T}|_{A}\to B]$, then $f\in[\mathcal{T}|_{A^{\prime}}\to B]$. In other words,

[TABLE]

If we combine two functions with provenance $A_{1}$ and $A_{2}$ respectively, then the resulting function has provenance $A_{1}\cup A_{2}$. This is shown in the following lemma.

Lemma 0.

Let $f_{1}\in[\mathcal{T}|_{A_{1}}\to B_{1}]$, $f_{2}\in[\mathcal{T}|_{A_{2}}\to B_{2}]$, $h\colon B_{1}\times B_{2}\to C$. Then

[TABLE]

This brings us to the main soundness result of Algorithm 1. Proposition 9 states that we can equip each CFG node $N$ with evaluation functions that depend on the addresses determined by prov in the sense of Definition 6. The evaluation functions directly map the trace $\mathtt{tr}$ to the values of variables at node $N$ such that the values agree with the small-step CFG semantics relation $\overset{\mathtt{tr}}{\longrightarrow}$. Importantly, while the operational semantics are defined for each individual trace, the evaluation functions work for all traces. In Fig. 1, concrete examples of these evaluation functions for a simple program can be found.

Proposition 0.

Let $G$ be the CFG of a program $S$ without while loop statements. For each node $N$ and variable $x$, there exists an evaluation function

[TABLE]

such that for all traces $\mathtt{tr}$, initial program state $\sigma_{0}$ as in Definition 1, and execution sequence

[TABLE]

we have

[TABLE]

Proof.

The CFG $G$ of a program without while loops is a directed acyclic graph and the statement is proven by induction on $G$. See Appendix B.2. ∎

Intuitively, $V^{N}_{x}$ computes the value of $x$ before executing $N$. These evaluation functions $V^{N}_{x}$ can be lifted to evaluation functions for expressions $V^{N}_{E}$:

[TABLE]

By Lemma 8, if $V^{N}_{y}\in[\mathcal{T}|_{\textnormal{prov}(N,y)}\to\mathcal{V}]$ for all $y\in\textnormal{vars}(E)$, then $V^{N}_{E}\in[\mathcal{T}|_{\textnormal{prov}(N,E)}\to\mathcal{V}]$.

3.3. Static Factorisation Theorem for Programs Without Loops

With Proposition 9, we can prove the first factorisation theorem by combining the evaluation functions $V^{N}_{x}$ to express the program density $p$ in a factorised form.

Theorem 10.

Let $G$ be the CFG for a program $S$ without while loop statements. Let $N_{1},\dots,N_{K}$ be all sample nodes in $G$. For each sample node ${N_{k}=\textnormal{Assign}(x_{k}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{k},f^{k}(E_{1}^{k},\dots,E_{n_{k}}^{k})))$, let

[TABLE]

Then, there exist functions $p_{k}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}]$ such that for all $\mathtt{tr}\in\mathcal{T}$, if $p_{S}(\mathtt{tr})\neq\textnormal{undefined}$, then

[TABLE]

Proof.

For each sample node $N_{k}$, define $b_{k}(\mathtt{tr})\coloneqq\bigwedge_{N^{\prime}\in\textnormal{BP}(N_{k})}t^{k}_{N^{\prime}}(V^{N^{\prime}}_{{\textnormal{condexp}(N^{\prime})}}(\mathtt{tr}))$ where $t^{k}_{N^{\prime}}(v)=v$ if $N_{k}$ is in the consequent branch of $N^{\prime}$ and $t^{k}_{N^{\prime}}(v)=\neg v$ if $N_{k}$ is in the alternate branch. By Lemma 8, we have $b_{k}\in[\mathcal{T}|_{\textstyle\bigcup_{N^{\prime}\in\textnormal{BP}(N_{k})}\textnormal{prov}(N^{\prime},\textnormal{condexp}(N^{\prime}))}\to\{\mathtt{true},\mathtt{false}\}].$

Here, $V^{N^{\prime}}_{\textnormal{condexp}(N^{\prime})}(\mathtt{tr})$ are precisely the evaluations of the branch conditions for node $N_{k}$. Thus, $b_{k}(\mathtt{tr})$ is $\mathtt{true}$ if $N_{k}$ is in the execution sequence for $\mathtt{tr}$ else $\mathtt{false}$. Define the factor $p_{k}$ as

[TABLE]

By Proposition 9, $V^{N_{k}}_{E_{i}^{k}}\in[\mathcal{T}|_{\textnormal{prov}(N_{k},E_{i}^{k})}\to\mathcal{V}]$ and $V^{N_{k}}_{E_{0}^{k}}\in[\mathcal{T}|_{\textnormal{prov}(N_{k},E_{0}^{k})}\to\textnormal{addresses}(N_{k})]$. Thus, it follows from Lemma 8 that $p_{k}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}]$ by construction. Lastly, from Proposition 9 we know that for the end-state $\sigma_{l}(\textbf{p})=V^{\text{END}}_{\textbf{p}}(\mathtt{tr})=p_{G}(\mathtt{tr})$. Since $b_{k}(\mathtt{tr})=\mathtt{true}$ if and only if $N_{k}$ is in the execution sequence for $\mathtt{tr}$ and since the $\textnormal{pdf}_{f^{k}}(\dots)$ term in (5) is precisely the value that gets multiplied to p at sample statements (1), we have $p_{G}(\mathtt{tr})=\prod_{k=1}^{K}p_{k}(\mathtt{tr})$. ∎

Again, consider the example program in Fig. 1. The density of this model factorises like below:

[TABLE]

This is exactly the factorisation obtained from Theorem 10. However, if we would replace the if-statement with a nonsensical one, (if (b == 1) then (m = 1) else (m = 1)), then the fourth factor is in $[\mathcal{T}|_{\{\textnormal{"x"},\textnormal{"b"},\textnormal{"s"}\}}\to\mathbb{R}_{\geq 0}]$ with $V^{N_{7}}_{\mathtt{m}}(\mathtt{tr})=\textnormal{ife}(\mathtt{tr}(\textnormal{"b"}),1,1)=1$, where we use ife as a short-hand for the if-else case distinction. In this case, the static analysis over-approximates the true provenance with the spurious dependency on ”b”.

3.3.1. Equivalence to Bayesian Networks

If we restrict the way sample statements are used in our PPL, we can establish a known equivalence to Bayesian networks. A Bayesian network (Koller and Friedman, 2009) is a directed acyclic graph $\mathcal{G}$ over a finite set of nodes $X_{1},\dots,X_{n}$ which represent random variables, such that the joint distribution factorises according to the graph $\mathcal{G}$,

[TABLE]

If we assume that all sample statements of a program have a constant unique address $\alpha_{k}$, ${N_{k}=\textnormal{Assign}(x_{k}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(\alpha_{k},f^{k}(E_{1}^{k},\dots,E_{n_{k}}^{k})))$, then we can identify each factor $p_{k}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}]$ of Theorem 10 with the unique address $\alpha_{k}\in\mathbf{Strings}$. We construct the Bayesian network $\mathcal{G}$ by mapping each address to a random variable $X_{\alpha_{k}}$. Since $\alpha_{k}$ is unique, $p_{k}$ is the only relevant factor for $X_{\alpha_{k}}$ and can be interpreted as conditional density function of $X_{\alpha_{k}}$. To see this, rewrite

[TABLE]

where the function $\mathbf{1}_{\mathtt{null}}(v)$ equals 1 if $v=\mathtt{null}$ else 0. If the node $N_{k}$ is not executed, then the value of $X_{\alpha_{k}}$ is assumed to be $\mathtt{null}$ and we add the $\mathbf{1}_{\mathtt{null}}\left(\mathtt{tr}(\alpha_{k})\right)$ term to $p_{k}$. As $\textnormal{pdf}_{f_{k}}$ is a probability density function and $\mathbf{1}_{\mathtt{null}}$ is the density of a Dirac distribution centered at $\mathtt{null}$, the factor $p_{k}$ is also a density function for each choice of values in $\mathtt{tr}$ at addresses $\alpha\in A_{k}\setminus\{\alpha_{k}\}$. Thus, we introduce the edge $\alpha_{j}\to\alpha_{k}$ if $\alpha_{j}\in A_{k}\setminus\{\alpha_{k}\}$, such that $\textnormal{parents}_{\mathcal{G}}(X_{\alpha_{k}})=A_{k}\setminus\{\alpha_{k}\}$ and

[TABLE]

In the above construction, it is important that each sample statement has a unique address. This guarantees no cyclic dependencies and a well-defined conditional probability distributions.

The Bayesian network equivalent to the example program of Fig. 1 is drawn below. Since the value of the program variable $\mathtt{m}$ depends on which program branch is taken during execution, the random variable $X_{\textnormal{"x"}}$ not only depends on $X_{\textnormal{"mu"}}$, but also on $X_{\textnormal{"b"}}$.

$\displaystyle P(X_{\textnormal{"b"}})$ $\displaystyle=\textnormal{pdf}_{\textnormal{Bernoulli}}(X_{\textnormal{"b"}};0.5)$

$\displaystyle P(X_{\textnormal{"s"}})$ $\displaystyle=\textnormal{pdf}_{\textnormal{InverseGamma}}(X_{\textnormal{"s"}};1,1)$

$\displaystyle P(X_{\textnormal{"mu"}}|X_{\textnormal{"b"}})$ $\displaystyle=\delta_{X_{\textnormal{"b"}}}\textnormal{pdf}_{\text{Normal}}(X_{\textnormal{"mu"}};0,1)+(1-\delta_{X_{\textnormal{"b"}}})\mathbf{1}_{\mathtt{null}}(X_{\textnormal{"mu"}})$

$\displaystyle P(X_{\textnormal{"x"}}|X_{\textnormal{"b"}},X_{\textnormal{"mu"}},X_{\textnormal{"s"}})$ $\displaystyle=\textnormal{pdf}_{\textnormal{Normal}}(X_{\textnormal{"x"}};\textnormal{ife}(X_{\textnormal{"b"}},X_{\textnormal{"mu"}},1),X_{\textnormal{"s"}})$

3.3.2. Equivalence to Markov Networks

If the sample addresses are not unique or not constant, then the factorisation is in general not a Bayesian network. However, we can show equivalence to Markov networks. A Markov network (Koller and Friedman, 2009) is an undirected graph $\mathcal{H}$ over random variables $\vv{X}$ that represents their dependencies. We consider Markov networks where the joint distribution of $\vv{X}$ factorises over a set of cliques $\mathcal{D}=\{D_{1},\dots,D_{n}\}\subseteq\textnormal{cliques}(\mathcal{H})$, with $\bigcup_{i=1}^{n}D_{i}=\vv{X}$:

[TABLE]

We again identify a random variable $X_{\alpha}$ for each address $\alpha\in\textstyle\bigcup_{k=1}^{K}A_{k}$ and construct a Markov network $\mathcal{H}$ by connecting node $X_{\alpha}$ to $X_{\beta}$ if $\alpha\in A_{k}\land\beta\in A_{k}$ for any $k$. Thus, $D_{k}=\{X_{\alpha}\colon\alpha\in A_{k}\}$ forms a clique and $p_{S}$ factorises over $\mathcal{H}$ with $\phi_{k}(D_{k})=p_{k}(\{\alpha\mapsto X_{\alpha}\colon\alpha\in A_{k}\})$.

Consider the program with constant but non-unique sample addresses in LABEL:lst:hurricane1. With nine sample statements we get a Markov network over $\{X_{\textnormal{"F"}},X_{\textnormal{"P0"}},X_{\textnormal{"P1"}},X_{\textnormal{"D0"}},X_{\textnormal{"D1"}}\}$. The nine factors correspond to the address sets $\{\textnormal{"F"}\}$, $\{\textnormal{"F"},\textnormal{"P0"}\}$, $\{\textnormal{"F"},\textnormal{"P1"}\}$, $\{\textnormal{"F"},\textnormal{"D0"},\textnormal{"P1"}\}$, $\{\textnormal{"F"},\textnormal{"D1"},\textnormal{"P0"}\}$,$2\times\{\textnormal{"F"},\textnormal{"P0"},\textnormal{"D0"}\}$, and $2\times\{\textnormal{"F"},\textnormal{"P1"},\textnormal{"D1"}\}$. In the consequent branch of the program, the distribution in the sample statement with address ”P1” depends on ”D0”, while in the alternative branch, the distribution of ”P0” depends on ”D1”. Thus, there is the dependency cycle $\textnormal{"P0"}\to\textnormal{"D0"}\to\textnormal{"P1"}\to\textnormal{"D1"}\to\textnormal{"P0"}$, which means the obtained factorisation is not a Bayesian network representation for the program. Note that the program can also be rewritten with dynamic addresses to eliminate the if statement. In this case, the factorisation still admits cyclic dependencies.

Lastly, while the programs considered up until this point were equivalent to finite Bayesian networks or finite Markov networks, note that the small program in LABEL:lst:poisson is equivalent to a Markov network with an infinite number of nodes $\{X_{\textnormal{"n"}}\}\cup\{X_{\textnormal{"x}\_i\textnormal{"}}\colon i\in\mathbb{N}\}$.

4. Static Factorisation for Programs With Loops

In this section, we describe how the argument for proving the factorisation theorem for programs without loops, can be extended to generalise the result to programs with loops.

4.1. Unrolled Control-Flow Graph

The translation rules for programs to their CFG can be extended to support while loops as follows. Let $G_{S}$ be the sub-CFG of statement $S$. The CFG of ${{(\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{while}}}}}}\;E\;\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{do}}}}}}\;S)$ is given below on the left:

  

For this CFG, $\textnormal{prov}(N,x)$ is still well-defined by Algorithm 1, but the proof of Proposition 9 does not work anymore since it requires a directed acyclic graph.

To establish Proposition 9 for programs with while loops, we introduce a second type of control-flow graph, the unrolled CFG, which again is a directed acyclic graph. All translation rules are as before except for while loops, which is given below and illustrated above on the right.

The resulting graph $G$ contains for every $i\in\mathbb{N}$ branch nodes $B_{i}=\textnormal{Branch}(E)$ and copies of the sub-CFG $G_{i}$ of statement $S$. It contains a single join node $J$ for which we have $(B_{i},J)\in\textnormal{BranchJoin}(G)$ for all $i$. Let $N_{\textnormal{first}}^{i}$ be successor of the start node in $G_{i}$ and $N_{\textnormal{last}}^{i}$ the predecessor of the end node in $G_{i}$. In the new graph $G$, the first branch node $B_{1}$ is successor of the start node and the end node is the successor of the join node $J$. The branch nodes $B_{i}$ have two successors $N_{\textnormal{first}}^{i}$ and $J$. $N_{\textnormal{last}}^{i}$, the last node of sub-graph $G_{i}$, resumes in the next branch node $B_{i+1}$.

This definition implies that the number of nodes of $G$ is countable. Furthermore, since the unrolled CFG consists only of start, end, branch, join, and assign nodes, we can equip it with the same semantics as the standard CFG. As before, the semantics are equivalent.

Proposition 0.

For all program $S$ with CFG $G$ and unrolled CFG $H$, it holds that for all traces $\mathtt{tr}$ and initial states $\sigma_{0}$

[TABLE]

Proof.

By structural induction. See Appendix B.3. ∎

As mentioned, Algorithm 1 still works for the unrolled CFG. However, since the unrolled CFG contains an infinite number of nodes, in practice, we run this algorithm on the standard CFG with a finite number of nodes avoiding loops by marking visited nodes. We will show that the computations on the standard CFG produce an over-approximation of the provenance set in the unrolled CFG. But first, we define the connection between the two graph types.

Definition 0.

Let $G$ be the CFG and $H$ the unrolled CFG of program $S$. For each node $M\in H$, we define $\iota_{\textnormal{cfg}}(M)\in G$ as the node from which $M$ was copied in the construction of $H$.

The definition of $\iota_{\textnormal{cfg}}$ is best understood when considering the example in Fig. 2, where the CFG and unrolled CFG of a program are depicted side-by-side and connected with the mapping $\iota_{\textnormal{cfg}}$.

The following lemma states that we can over-approximate the provenance set for variable $x$ at node $M$ in the unrolled CFG $H$ by applying Algorithm 1 to the corresponding node $\iota_{\textnormal{cfg}}(M)$ in $G$.

Lemma 0.

Let $G$ be the CFG and $H$ the unrolled CFG of program $S$. For each node $M\in H$ following statements hold:

[TABLE]

[TABLE]

[TABLE]

Proof.

See Appendix B.4. ∎

4.2. Static Factorisation Theorem for Programs With Loops

The unrolled CFG has the important property that it is a directed acyclic graph with a potentially infinite number of nodes, but with a single root node. Further, each node $N\neq\textnormal{START}$ still has exactly one predecessor except for join nodes. Therefore, Proposition 9 can be generalised to unrolled CFGs with only a slight modification to the proof.

Proposition 0.

Let $H$ be the unrolled CFG of a program $S$. For each node $M$ and variable $x$, there exists an evaluation function

[TABLE]

such that for all traces $\mathtt{tr}$, initial program state $\sigma_{0}$ as in Definition 1, and execution sequence

[TABLE]

we have

[TABLE]

Proof.

The unrolled CFG $H$ is constructed in a way such that we can prove the statement by well-founded induction on the relation $M_{1}\prec M_{2}\Leftrightarrow M_{1}\textit{ is parent of }M_{2}$ (i.e. there is an edge from $M_{1}$ to $M_{2}$ in $H$) similar to the proof of Proposition 9. The full proof is given in Appendix B.5. ∎

As before, Proposition 4 allows us to explicitly construct the factorisation of the program density.

Theorem 5.

Let $G$ be the CFG for a program $S$. Let $N_{1},\dots,N_{K}$ be all sample nodes in $G$. For each sample node ${N_{k}=\textnormal{Assign}(x_{k}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{k},f^{k}(E_{1}^{k},\dots,E_{n_{k}}^{k})))$, let

[TABLE]

Then, there exist functions $p_{k}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}]$ such that for all $\mathtt{tr}\in\mathcal{T}$, if $p_{S}(\mathtt{tr})\neq\textnormal{undefined}$, then

[TABLE]

Proof.

We give a proof sketch and refer to Appendix B.6 for the full proof. For each sample node $M_{j}$ in the unrolled CFG $H$, we construct $\tilde{p}_{j}\in[\mathcal{T}|_{\tilde{A}_{j}}\to\mathbb{R}_{\geq 0}]$ as in the proof of Theorem 10. Next, we group the nodes $M_{j}$ by their corresponding CFG node: $p_{k}=\prod_{j\colon\iota_{\textnormal{cfg}}(M_{j})=N_{k}}\tilde{p}_{j}$. From Lemma 3, we have that if $\iota_{\textnormal{cfg}}(M_{j})=N_{k}$, then $\tilde{A}_{j}\subseteq A_{k}$. Thus, $p_{k}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}]$. Like in the proof of Theorem 10 it follows from Proposition 4 that $p_{S}(\mathtt{tr})=p_{G}(\mathtt{tr})=\prod_{k=1}^{K}p_{k}(\mathtt{tr})$. ∎

4.2.1. Equivalence to Markov Networks

As in Section 3.3.2, Theorem 5 implies that a program $S$ with $K$ sample statements is equivalent to a Markov network $\mathcal{H}$ with nodes $\{X_{\alpha}\colon\alpha\in\textstyle\bigcup_{k=1}^{K}A_{k}\}$ even if it contains loops. Again, the sets $D_{k}=\{X_{\alpha}\colon\alpha\in A_{k}\}$ form cliques in $\mathcal{H}$ and $p_{S}$ factorises over $\mathcal{H}$ with $\phi_{k}(D_{k})=p_{k}(\{\alpha\mapsto X_{\alpha}\colon\alpha\in A_{k}\})$: $P(\vv{X})=\prod_{k=1}^{K}\phi_{k}(D_{k}).$

5. Optimising Posterior Inference Algorithms

While the statically obtained factorisation is provably correct, it is natural to ask whether the quality of the factorisation is good enough to be of practical use. To answer this question, we have developed a program transformation based on the static factorisation that generates sub-programs for each factor (described in Section 5.1). This program transformation operates on probabilistic programs that are implemented in our research PPL which models the formal language introduced in Section 2. It is implemented in Julia and can be viewed as a subset of PPLs like Gen, Turing, or Pyro. Many posterior inference algorithms can be improved by exploiting the factorisation of the model density. In Section 5.2, we demonstrate that if the change to the current program trace is small as in the light-weight Metropolis Hastings kernel (Wingate et al., 2011), then we can leverage the sub-programs to compute the acceptance probability faster. In Section 5.3, we show that we can reduce the variance of the ELBO gradient estimation in variational inference. Lastly, in Section 5.4, we compile sub-programs to enable an iterative implementation of sequential Monte Carlo. The realised optimisation tricks are well-known but typically restricted to programs with a finite number of random variables, PPLs that require the user to explicitly construct the graphical structure, or dynamic approaches. In contrast, our novel transformation operates on the source code and makes the optimisations statically available for programs containing while loops and dynamic addresses. We evaluate our implementation (Böck, 2025) on a benchmark set of 16 probabilistic programs, the majority of which make use of while loops to declare an unbounded number of random variables.

5.1. Generating Sub-programs from the Static Factorisation

We generate a sub-program for each sample statement (and thus for each factor) by deriving a program slicing method from the static provenance analysis. We want to generate a sub-program for each sample node ${N_{k}=\textnormal{Assign}(x_{k}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{k},f^{k}(E_{1}^{k},\dots,E_{n_{k}}^{k})))$ that

1. (1)

enables us to continue execution directly from $N_{k}$ from any state $\sigma$, 2. (2)

stops execution once all dependencies of $N_{k}$ were reached, 3. (3)

supports multiple inference algorithms.

To realise (1) and (2), we collect all sample nodes that depend on $N_{k}$. From Theorem 5 we know that a sample node $N_{j}$ depends on $N_{k}$ if $N_{k}$ contributes to the provenance $A_{j}$ of the factor $p_{j}\in[\mathcal{T}|_{A_{j}}\to\mathbb{R}_{\geq 0}]$ as defined in Eq. 9. As explained in Section 3.2, these nodes can be found by modifying Algorithm 1 on line 8 to collect sample nodes instead of their addresses. For this modified algorithm provnode, we have $\textnormal{prov}(N,x)=\bigcup_{N^{\prime}\in\textnormal{provnode}(N,x)}\textnormal{addresses}(N^{\prime})$ and $N_{j}$ depends on $N_{k}$ if any variable $y$ in any expression $E_{i}^{j}$ in the sample statement of $N_{j}$ or in the conditional expression of any branch parent of $N_{j}$ depends on $N_{k}$. We write $j\in\mathcal{J}_{k}$ if $N_{j}$ depends on $N_{k}$.

We slice the program to only retain nodes that lie on a path from $N_{k}$ to $N_{j}$ for some $j\in\mathcal{J}_{k}$. More formally, let $G_{k}$ be the CFG that contains node $N$ if there exists $j\in\mathcal{J}_{k}$ such that there is a path from $N_{k}$ to $N_{j}$ that visits $N$ and does not revisit $N_{k}$. The start node of $G_{k}$ has successor $N_{k}$ and the nodes $N_{j}$, $j\in\mathcal{J}_{k}$, are connected to the end node of $G_{k}$ if they have no successor in $G_{k}$. Lastly, sample nodes $N_{j}$ in $G_{k}$ that do not depend on $N_{k}$, $j\notin\mathcal{J}_{k}$, are transformed to read nodes for which ${(\sigma,x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{read}}}}}}(E_{0}))\Downarrow^{\mathtt{tr}}\sigma[x\mapsto\mathtt{tr}(\sigma(E_{0}))]$. This construction is made precise in Appendix D.1, where we also prove following correctness theorems as consequences of Theorem 5.

Theorem 1.

Let $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ be traces that only differ at address $\alpha$, i.e. $\mathtt{tr}_{1}(\alpha)\neq\mathtt{tr}_{2}(\alpha)$ and $\forall\beta\neq\alpha\colon\mathtt{tr}_{1}(\beta)=\mathtt{tr}_{2}(\beta).$ Let the execution sequences in the unrolled CFG $H$ be equal up to sample node ${M_{\alpha}=\textnormal{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},f(E_{1},\dots,E_{n})))$ and state $\sigma$, such that $\alpha=\sigma(E_{0})$:

[TABLE]

If $M_{\alpha}$ is the only sample node in the execution sequences whose address expression evaluates to $\alpha$, then for the sliced CFG $G_{k}$ of sample node $N_{k}=\iota_{\textnormal{cfg}}(M_{\alpha})$, it holds that

[TABLE]

Proof.

Follows from Theorem 5 and two technical lemmas. See Appendix D.2. ∎

Corollary 0.

Assume that for each trace $\mathtt{tr}$ there is at most one sample node $M_{\alpha}$ in its execution sequences such that the address expression of $M_{\alpha}$ evaluates to $\alpha$. For initial program state $\sigma_{0}$ as in Definition 1, define

[TABLE]

Then, $\Delta_{\alpha}(\mathtt{tr})\coloneqq\log p_{G}(\mathtt{tr})-\log p_{\alpha}(\mathtt{tr})$ is independent of address $\alpha$, $\Delta_{\alpha}\in[\mathcal{T}|_{\mathbf{Strings}\setminus\{\alpha\}}\to\mathbb{R}_{\geq 0}].$

Proof.

Follows directly from Theorem 1. See Appendix D.3. ∎

Lastly, to fulfill (3), we further transform the program to take as input an execution context ctx. The context controls the behavior when executing a sample statement, e.g. to resample a variable or to accumulate the density. To actually handle the effects of different inference algorithms, we replace sample statements in the original program with three abstract methods that delegate the behavior to the context ctx. In $G_{k}$, the statement of $N_{k}$ is transformed to visit, the statements of the nodes that depend on $N_{k}$ to score, and all other sample nodes that are on a path between $N_{k}$ and some $N_{j}$, $j\in\mathcal{J}_{k}$, but do not depend on $N_{k}$, are transformed to read. Concrete instances of such contexts will be presented in the following sections and formalised in Appendix D.

We have implemented the sub-program generation as a compilation to the Julia programming language making use of multiple dispatch for contextualised execution. To enable continuation from any program state $\sigma$ we also inject a mutable state struct $\rho$. As the CFG slicing implementation is a fairly standard graph traversal algorithm, we refer to our implementation (Böck, 2025) for further details like the translation from CFGs to source code.

The sub-program generation process is best illustrated with three examples. Consider program 1 in Figure 3. To generate a sub-program for the factor of sample statement ”B”, we omit the sample statements ”A” and ”E” as they do not lie on a path from ”B” to ”D” — the only sample statement that depends on ”B”. However, we still need to read the value of program variable A from the cached program state $\rho$ and we read the value for random variable ”C” from the trace. Continuing execution directly from a sample statement by caching the program state means that we can continue a loop from a specific iteration as exemplified by program 2 in Figure 3. This is particularly powerful, if the only dependency of a sample statement belongs to the same loop iteration. In this case, an entire while loop can be replaced by a single iteration which is shown for program 3 in Figure 3. Note that this form of loop-iteration independence is not captured by the factorisation of Proposition 4. Instead, it is product of the program slicing implementation.

5.2. Speeding-up Acceptance Probability Computation for Single-Site Trace Updates

In this section we introduce execution contexts to speed-up the single-site or light-weight Metropolis Hastings inference (LMH) algorithm (Wingate et al., 2011). The LMH algorithm works by selecting a single address $\alpha$ in each iteration at random and proposing a new value $v$ for the current trace $\mathtt{tr}(\alpha)$. Even though we select a single address for the update, the proposal may have to produce additional values for sample statements that were not executed for trace $\mathtt{tr}$, but due to stochastic branching or dynamic addressing are now executed for the new trace $\mathtt{tr}^{\prime}$. For instance in Figure 3 program 2, if we change the value of the last random variable $\textnormal{"b}\_i\textnormal{"}$ from 0 to 1, then we may have to sample values for multiple $\textnormal{"b}\_j\textnormal{"}$ where $j>i$ until we get $\mathtt{tr}^{\prime}(\textnormal{"b}\_j\textnormal{"})=0$.

The proposed trace $\mathtt{tr}^{\prime}$ is accepted (made current) with probability $A=\min\left(1,\frac{p(\mathtt{tr}^{\prime})\times Q_{\alpha}(\mathtt{tr}|\mathtt{tr}^{\prime})}{p(\mathtt{tr})\times Q_{\alpha}(\mathtt{tr}^{\prime}|\mathtt{tr})}\right),$ where $Q_{\alpha}$ captures the probability of the proposal process (sampling a value for $\alpha$ and all additional values as explained above). Then, the set of all generated traces is an approximation to the posterior distribution of the model.

In our approach, if the address comes from sample node $N_{k}$, then we can optimise the runtime of this algorithm by noting that $\log p(\mathtt{tr}^{\prime})-\log p(\mathtt{tr})=\sum_{j=1}^{K}\log p_{j}(\mathtt{tr}^{\prime})-\log p_{j}(\mathtt{tr})=\sum_{j\in\mathcal{J}_{k}\cup\{k\}}\log p_{j}(\mathtt{tr}^{\prime})-\log p_{j}(\mathtt{tr})$. Thus, by being able to compute the factors $p_{j}(\mathtt{tr}^{\prime})$ and $p_{j}(\mathtt{tr})$ separately, we can compute the ratio of densities $p(\mathtt{tr}^{\prime})/p(\mathtt{tr})$ by summing fewer terms.

By designing execution contexts for the sub-programs to handle the effect of the LMH kernel, we statically enable this optimisation for programs with an unbounded number of random variables and dynamic dependencies. In the following we informally describe the implementation of two execution contexts for resampling address $\alpha$ at sample node $N_{k}$. A formal treatment can be found in Appendix D.1.

When executed with ForwardCtx, the sub-program is evaluated with respect to the proposed trace $\mathtt{tr}^{\prime}$. At visit(ForwardCtx,$\rho$,$E_{0}$,$f(\dots$)) we sample a new value $v$ for address $\alpha=\rho(E_{0})$ from some proposal distribution $q$ and return $v$. As before, $\rho(E_{0})$ denotes $E_{0}$ evaluated with respect to $\rho$. In the background, we accumulate the log-probability $\Delta P\leftarrow\Delta P+\log\textnormal{pdf}_{f}(v;\dots)$ and proposal log-probability $\Delta Q\leftarrow\Delta Q-\log\textnormal{pdf}_{q}(v;\dots)$. At the start, $\Delta P$ and $\Delta Q$ are set to [math]. At score(ForwardCtx,$\rho$,$E_{0}$,$f(\dots$)) we propose a new value $v$ only if the address $\rho(E_{0})$ is not sampled for the current trace. We always update $\Delta P$ as before and only update $\Delta Q$ if we have proposed a new value. At read(ForwardCtx,$\rho$,$E_{0}$), we simply return $\mathtt{tr}(\rho(E_{0}))$. At each of the above statements we cache the program state in $\rho$. This introduces book-keeping cost, but as we will see in Section 5.2.1, the cost is outweighed by the performance gained from only running sub-programs.

When executed with BackwardCtx, the sub-program is evaluated with respected to the current trace $\mathtt{tr}$. At visit(BackwardCtx,$\rho$,$E_{0}$,$f(\dots$)) we update $\Delta P\leftarrow\Delta P-\log\textnormal{pdf}_{f}(v;\dots)$ and proposal log-probability $\Delta Q\leftarrow\Delta Q+\log\textnormal{pdf}_{q}(v;\dots)$ with the current value $v=\mathtt{tr}(\rho(E_{0}))$. At score(BackwardCtx,$\rho$,$E_{0}$,$f(\dots$)) the same $\Delta P$-updates are always performed while the $\Delta Q$-updates are only performed if the address $\rho(E_{0})$ is not sampled for the proposed trace. At statement read(BackwardCtx,$\rho$,$E_{0}$) we again return $\mathtt{tr}(\rho(E_{0}))$. The program state does not have to be cached when executing with BackwardCtx.

Lastly, after executing the sub-program in both contexts, we have $A=\min(1,\exp(\Delta P+\Delta Q))$. We accept the proposed trace with probability $A$. In this case, we perform some book-keeping to update the inference state to the proposed trace with the cached program states. Building on the soundness of our provenance analysis and Theorem 1, we proof correctness of this approach in Appendix D.4:

Corollary 0.

Assume that the address expression of at most one sample node in the execution sequence evaluates to $\alpha$ for all traces. Let $\mathtt{tr}^{\prime}$ be the LMH proposal at address $\alpha$ for $\mathtt{tr}$. Then, the execution sequences of $\mathtt{tr}$ and $\mathtt{tr}^{\prime}$ are equal up to $(\sigma,M_{\alpha})$, where $M_{\alpha}$ is the first and only sample node whose address expression evaluates to $\alpha$. For the sliced CFG $G_{k}$, where $\iota_{\textnormal{cfg}}(M_{\alpha})=N_{k}$, we have

[TABLE]

The quantities $\log\bar{p}_{G_{k}}(\mathtt{tr}^{\prime},\sigma)$, and $\log Q_{\alpha}(\mathtt{tr}^{\prime}|\mathtt{tr})$ are computed by ForwardCtx and the quantities $\log\bar{p}_{G_{k}}(\mathtt{tr},\sigma)$ and $\log Q_{\alpha}(\mathtt{tr}|\mathtt{tr}^{\prime})$ are computed by BackwardCtx. The new trace $\mathtt{tr}^{\prime}$ can be constructed from $\mathtt{tr}$ and the output of both contexts.

5.2.1. Experiment Results

To evaluate the speed-up made possible by our static factorisation, we gathered a benchmark set of 16 probabilistic programs. We specifically collected probabilistic programs from prior research that make use of while loops to define an unbounded number of random variables. References for the programs can be found in our replication package (Böck, 2025). Only 5 out of the 16 models define just a finite number of random variables with fixed dependency structure. It takes less than five seconds to generate the sub-programs for all models.

For each model, we run two implementations of the LMH inference algorithm. As a baseline, the standard implementation of LMH computes the acceptance probability $A$ by computing the density $p(\mathtt{tr}^{\prime})$ from scratch. Our optimised version of LMH leverages the sub-programs to compute $\log p(\mathtt{tr}^{\prime})-\log p(\mathtt{tr})$ when calculating $A$. For the five finite models with fixed structure, we additionally implemented a factorisation that is representative of PPLs that explicitly construct a probabilistic graphical model. This factorisation enables LMH to exploit the full dependency structure of the model and serves as comparison. We have tested correctness of our LMH implementations by asserting that exactly the same acceptance probabilities and samples are computed.

In Table 1, we report the results of the experiment run a M2 Pro CPU. For each model and implementation, we ran LMH with 10 random starting traces for a total of $N$ samples and average the measurements over 10 experiment repetitions.

The first observation that we want to highlight is that computing $\log p(\mathtt{tr}^{\prime})-\log p(\mathtt{tr})$ with our approach can lead to longer runtime if there is no structure to exploit in the model. Consider the linear regression model in which the 100 observations depend on the slope and intercept variables. When changing the value of one of these variables, computing $\log p(\mathtt{tr}^{\prime})-\log p(\mathtt{tr})$ requires 100 additions and 100 subtractions, while $\log p(\mathtt{tr}^{\prime})$ only takes 100 additions. The same is true for the Marsaglia model which implements a rejection sampling method for the Normal distribution such that all variables depend on each other.

The factorisation approach shines when there is a lot of structure in the models like for the Aircraft, Captcha, Dirichlet Process, and Urn models, where the runtime is decreased up to a factor of 2. It is particularly powerful, when there are latent variables that only affect independent data variables. This is the case for the Gaussian mixture model (GMM) where there is one cluster allocation variable per data point and for the Latent Dirichlet Allocation model (LDA) where there is similar dependency structure. For these models we observe 8-10x speed-ups. Note that the handwritten ”fixed-finite” factorisations of the models are faster than the ”static analysis” factorisations, because LMH has to do less book-keeping for them. They do not leverage additional dependency structure.

The only model where our approach fails to effectively leverage the true dependency structure is the Hidden Markov Model (HMM). Even though the current state only affects the next state, our static analysis over-approximates this structure and outputs that the current state affects all future states. Now, if we modify the program by unrolling the while loop, i.e. repeating its body fifty times, our factorisation approach results in an improved runtime of 2.605 $\mu s$ (13x speed-up), because in this case the dependency structure is not over-approximated anymore due to loops. Note that the PCFG model essentially has the same structure as HMM, but can transition to a terminal state at every loop iteration. In this case, indeed the current state affects all future states and our static analysis is exact, but no significant improvements in runtime are measurable.

Lastly, we note that the Hurricane model is the only model in our benchmark with a finite number of random variables but a dynamic dependency structure (see Section 3.3.2). This model can also not be represented in a conventional Bayesian network structure, but our method is applicable and can slightly reduce the runtime.

5.2.2. Comparison to Related Methods

The approach presented in Section 5 is the first to generate optimised MCMC kernels of models with unbounded number of random variables via a static program transformation. Importantly, its correctness follows from the theoretical foundation established in Section 4. Nevertheless, there exist methods that aim to accelerate MCMC kernels for a similar class of models with other approaches. In this section, our goal is to make a fair comparison to those methods by focusing on the relative speed-ups achieved. A direct runtime comparison is not meaningful due to differences in programming languages and the implementation of inference algorithms. More related approaches are described in Section 6.

First, Ritchie et al. (2016) presented the C3 approach to enable incrementalised re-execution of Metropolis Hastings kernels using continuations and function call caching. To be effective this methods requires programs to be implemented with recursion instead of while loops. We re-implemented the models in recursive form in WebPPL (Goodman and Stuhlmüller, 2014) to report the relative speed-up achieved with C3 in Table 2. Our approach achieves better speed-ups for all models except for the HMM. As explained before, our static analysis greatly over-approximates the dependency structure of this model. Note that Ritchie et al. (2016) report even more significant speed-ups for the GMM, LDA, and HMM models for larger data-set sizes. Our approach scales the same for GMM and LDA.

Further, modern probabilistic programming systems like Gen (Cusumano-Towner et al., 2019) include language constructs like Map, Unfold, and Recurse to model independence structure. Thus, we re-implemented LMH and the 16 models only in Gen with the so-called combinator constructs and report the results in Table 2. Gen’s inference backend can only leverage the combinator constructs if the root model function is written in the static modelling language – a non-universal subset of the language with no control statements. This effectively abstracts away loops as Map, Unfold, or Recurse nodes in an internal graphical representation. Note that at the time of writing Gen does not implement incremental computation for the Recurse combinator, a construct required for models where the condition of the while loop depends on the body like the Geometric model. We cannot apply our general LMH implementation to those models and mark them with ”-”. Overall, the order of magnitude of speed-up achieved with combinators is similar to our static approach. There seems to be more overhead in Gen as it supports general programmable inference in contrast to our research PPL. We emphasise that our approach achieves these speed-ups fully automatically from source code without the need to explicitly model independence structure with language constructs.

5.3. Reducing the Variance of an ELBO Gradient Estimator in Variational Inference

Variational inference takes a different approach to MCMC methods by fitting a so-called variational distribution $Q_{\phi}$ to the posterior, where $\phi$ is a set of parameters. The objective is to maximise the evidence lower-bound $\textnormal{ELBO}(\phi)=\mathbb{E}_{z\sim Q_{\phi}}{\left[\log P(z)-\log Q_{\phi}(z)\right]}$ via gradient ascent, where $P$ denotes a potentially unnormalised density. In this section, we consider Black-Box Variational Inference (BBVI) (Ranganath et al., 2014), which estimates the gradient of the ELBO with following identity $\nabla_{\phi}\textnormal{ELBO}(\phi)=\mathbb{E}_{z\sim Q_{\phi}}{\left[\nabla_{\phi}\log Q_{\phi}(z)\times(\log P(z)-\log Q_{\phi}(z))\right]}$. Unfortunately, this estimator often exhibits high variance. To reduce the variance Ranganath et al. (2014) use a mean-field variational distribution $Q_{\phi}(z)=\prod_{i=1}^{n}q_{i}(z_{i}|\phi_{i})$ to derive the improved estimator

[TABLE]

Note that the gradient is only with respect to $\phi_{i}$ and the expectation is with respect to the Markov blanket $z_{(i)}$ of $z_{i}$. Including $z_{i}$, this is the set of variables that directly influence $z_{i}$ according to $P$, variables that are directly influenced by $z_{i}$ – children of $z_{i}$ – and variables that directly influence a child of $z_{i}$. Then, $P^{i}$ is the density factoring out all other variables and $Q^{i}_{\phi}$ the variational distribution restricted to $z_{(i)}$. This reduces the variance and we will replicate this trick with our factorisation.

In our sub-program setting, we achieve the same optimisation by defining the execution context BBVICtx, which we formalise in Appendix D.5. First, assume that a trace $\mathtt{tr}$ is generated from a variational distribution, where each address $\alpha$ is sampled from $\mathtt{tr}(\alpha)\sim q_{\alpha}(\mathtt{tr}(\alpha)|\phi_{\alpha})$. The BBVICtx executes the program with respect to $\mathtt{tr}$. Similar to the ForwardCtx, it accumulates $\Delta\leftarrow\Delta+\log\textnormal{pdf}_{f}(v;\dots)$ at score(BBVICtx,$\rho$,$E_{0}$,$f(\dots$)) and returns $v=\mathtt{tr}(\rho(E_{0}))$ at read(BBVICtx,$\rho$,$E_{0}$)). At visit statement visit(BBVICtx,$\rho$,$E_{0}$,$f(\dots$)) we additional substract the proposal density $\Delta\leftarrow\Delta+\log\textnormal{pdf}_{f}(v;\dots)-\log q_{\alpha}(v|\phi_{\alpha})$, where we have $\alpha=\rho(E_{0})$ (compare Eq. 10). In Appendix D.5. we prove following correctness result as a consequence of Corollary 2:

Corollary 0.

Let $Q_{\phi}$ be a variational distribution over traces as formalised in Appendix D.5. Assume that the address expression of at most one sample node in the execution sequence evaluates to $\alpha$ for all traces. Under measurability assumptions given in Appendix D.5 and with $p_{\alpha}$ as in Corollary 2,

[TABLE]

The quantity $\log p_{\alpha}(\mathtt{tr})-\log q_{\alpha}(\mathtt{tr}(\alpha)|\phi_{\alpha})$ is computed by BBVICtx.

Note that in contrast to Equation 10, where the expectation is taken explicitly with respect to the Markov blanket $z_{(i)}$, in Corollary 4 the dependencies are accounted for via $p_{\alpha}(\mathtt{tr})$ (computed with a sub-program). Replacing $p$ in the integrand with $p_{\alpha}$ reduces the variance, because $p_{\alpha}$ is constant with respect to addresses that are independent from a change to the trace at address $\alpha$.

5.3.1. Experiment Results

To evaluate the approach we have implemented our BBVI optimisation and compare the improved estimator against the standard estimator with respect to the average gradient variance across all parameters. In 10 experiment repetitions, we performed 1000 optimisation steps and used 100 samples to estimate the gradient in each step. In Table 3, we report our results, where we omit programs with random variables that have stochastic support dimension and are incompatible with BBVI.

We compare against Pyro’s TraceGraph_ELBO variance reduction technique, which follows Schulman et al. (2015) and makes use of dynamic data provenance analysis. This comes with the benefit of having no over-approximation issues, but also comes with the drawback of not being able to correctly account for control flow dependencies. For Pyro we decreased the number of samples to 100 and number of repetitions to 1 due to long runtime. In Table 3, we see that for programs with only data dependencies we achieve comparable variance reduction with lower runtime penalty. Furthermore, with our approach we were able to reduce the variance for two models with control flow dependencies Aircraft and Marsaglia. Note that for the same reasons as in LMH, the static optimisation is only effective for the unrolled version of the HMM.

5.4. Iteratively Running Sequential Monte Carlo

Lastly, we briefly describe how our slicing technique can be adapted to the setting of sequential Monte Carlo (SMC)  (Chopin et al., 2020), which yet again takes a different approach to inference. In SMC a number of traces $\mathtt{tr}_{n}$, so-called particles, are evaluated in parallel. Data points are added iteratively to the inference process and the particles are resampled according to their data likelihood. More precisely at time $t$ we have data set $\mathcal{D}_{1:t}$ and define the corresponding density $p_{t}$. At inference step $t$, the particle $\mathtt{tr}$ is logarithmically weighted according to $\log p_{t}({\mathtt{tr}_{n}})-\log p_{t-1}({\mathtt{tr}_{n}})$ and the particles are resampled (duplicated or discarded) with respect to these weights.

A ”naive” implementation of SMC would compute $\log p_{t}({\mathtt{tr}_{n}})$ from scratch, while an iterative implementation is able to compute $\log p_{t}({\mathtt{tr}_{n}})-\log p_{t-1}({\mathtt{tr}_{n}})$ directly. In probabilistic programming this is often realised by some mechanism that enables the interruption and continuation of program execution. We enable this optimisation by slightly adapting the slicing method of Section 5.1 to generate sub-programs which compute the difference directly.

Namely, for sample node $N_{k}$ let $G_{k}^{\texttt{SMC}}$ be the CFG that contains node $N$ if there exists a sample node $N_{j}$ such that $N$ lies on a path from $N_{k}$ to $N_{j}$ which contains no other sample nodes. Again, $N_{k}$ will be the successor of the start node of $G_{k}^{\texttt{SMC}}$, but all other samples nodes will have the end node of $G_{k}^{\texttt{SMC}}$ as sole successor. As we only care about continuing execution from $N_{k}$ to the next sample node, this construction is independent of the provenance analysis and Theorem 5.

In our implementation, we again make use of an evaluation context which simply accumulates the density like the standard semantics. In Appendix D.6, we prove correctness:

Proposition 0.

Let $G$ be a CFG such that for initial program state $\sigma_{0}$ and trace $\mathtt{tr}$, we have execution sequence $(\sigma_{0},\text{START})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma_{i},M_{i})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma_{l},\text{END})$ in its unrolled version $H$. Let $(\sigma^{\prime}_{1},M^{\prime}_{1})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{2},M^{\prime}_{2})$ be a sub-sequence such that $\iota_{\textnormal{cfg}}(M^{\prime}_{1})=N_{k}$ and $\iota_{\textnormal{cfg}}(M^{\prime}_{2})=N_{j}$ are the only sample nodes. Then, the execution sequence in the unrolled version of the sliced CFG $G_{k}^{\texttt{SMC}}$ matches the sub-sequence: $(\sigma^{\prime}_{1},\text{START})\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{1},M^{\prime}_{1})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{2},M^{\prime}_{2})\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{2},\text{END}).$

Proposition 5 allows us to evolve the execution sequence of a particle $\mathtt{tr}_{n}$ from one sample node $N_{t}^{n}=N_{k}$ to the next $N_{t+1}^{n}=N_{j}$ by executing the corresponding sliced CFG $G_{k}^{\texttt{SMC}}$. The weight of the particle is updated with $\log p_{t}(\mathtt{tr}_{n})-\log p_{t-1}(\mathtt{tr}_{n})=\log\sigma_{2}^{\prime}(\textbf{p})-\log\sigma_{1}^{\prime}(\textbf{p})$ and particles may be resampled even if the execution sequences were stopped at different sample nodes $N_{t}^{n}$ (Lundén et al., 2021).

5.4.1. Experiment Results

To measure the performance gains achievable with our static approach, we compare the runtime of the iterative implementation of SMC based on the generated sub-programs against a naive implementation which computes $\log p_{t}(\mathtt{tr})$ from scratch. In Table 4 we report the speed-ups on the subset of models which can make use of iteratively adding data points in inference. For all models we set the number of particles to 100 and repeated the experiment 10 times. The measured speed-ups are significant and range from 3x - 10x. For reference, we benchmarked the runtimes achieved with WebPPL’s continuation passing style (CPS) approach against a naive SMC implementation in JavaScript. For most models this resulted in greater speed-ups compared to our static approach. This is due to the fact that, like in Section 5.2, in our approach computational overhead is introduced by caching and copying programs states when resampling particles.

6. Related Work

6.1. Semantics of Probabilistic Programming Languages

One of the first works to formally define the semantics of a probabilistic programming language was Kozen (1979). Derivatives of this semantics can be found throughout literature, for example see Hur et al. (2015), Dahlqvist and Kozen (2019), or Barthe et al. (2020). We give a brief (and incomplete) overview of existing semantic models for PPLs.

In recent years, a lambda calculus for probabilistic programming was presented by Borgström et al. (2016). This work and a domain theoretic approach for probabilistic programming (Vákár et al., 2019) had an influence on the development of the statistical PCF (Mak et al., 2021) (programming computable functions), which was used to prove differentiability of terminating probabilistic programs. Probabilistic programming was also examined from a categorical perspective (Heunen et al., 2017) and the semantics of higher order probabilistic programs is also researched (Dahlqvist and Kozen, 2019; Staton et al., 2016).

As mentioned, the semantics presented in this work generalise those of Gorinova et al. (2019) which formalises the core language constructs of Stan by interpreting a program as a function over a fixed set of finite variables. In our semantics, we interpret programs as functions over traces which are mappings from addresses to values similar to dictionaries. Another dictionary based approach to probabilistic programming can be found in Cusumano-Towner et al. (2020) and denotational semantics based on trace types are introduced in Lew et al. (2019). To the best of our knowledge there are only two existing semantic models that support user-defined addresses (Lew et al., 2023; Lee et al., 2019). Note that denotational semantics of a simplified version of the address-based PPL Gen are given in Cusumano-Towner (2020).

Furthermore, the connection of simple probabilistic programs to Bayesian networks is known (van de Meent et al., 2018). Borgström et al. (2011) consider a small imperative language and give semantics in terms of factor graphs – a concept related to Bayesian networks. Sampson et al. (2014) compile probabilistic programs to Bayesian networks by treating loops as black-box functions in order to verify probabilistic assertion statements. Previous work extended Bayesian networks to so-called contingent Bayesian networks in order to support models with an infinite number of random variables or cyclic dependencies (Milch et al., 2005), but they cannot be straightforwardly applied to PPLs with user-labelled sample statements and while loops. Paquet (2021) presented a interpretation of probabilistic programs in game semantics, which can be seen as a refinement of Bayesian networks. Most recently, Faggian et al. (2024) introduced a higher-order functional language for discrete Bayesian networks.

6.2. Graphical Representations, Optimisations, and Dependency Analysis

The advantages of representing a probabilistic model graphically are well-understood (Koller and Friedman, 2009). For this reason, some PPLs like PyMC (Salvatier et al., 2016), Factorie (McCallum et al., 2009), or Infer.NET (Minka, 2012) require the user to explicitly define their model as a graph. This is not possible if we want to support while loops in our programs. There are also other graphical representations that allow for optimised inference like binary decision diagrams (Holtzen et al., 2020) or sum-product networks (Saad et al., 2021). Venture (Mansinghka et al., 2018) represents traces graphically to facilitate efficient recomputation of subparts of the program. Similar to Gen (Cusumano-Towner et al., 2019), Pyro (Bingham et al., 2019) implements the plate and markov constructs to model independence and reduce variance in variational inference.

Modifying the Metropolis Hastings algorithm to exploit the structure of the model is also common practice. Wu et al. (2016) optimised BLOG to dynamically keep track of the Markov blanket of a variable among other improvements. For their Metropolis Hastings implementation, they reported significant speed-ups for the Hurricane model, instances of the Urn model, and GMMs with unknown number of clusters. The Shred system (Yang et al., 2014) traces Church (Goodman et al., 2012) programs and slices the them to identify the minimal computation needed to compute the MH acceptance probability, and translates them to C++. However, their approach comes with high compilation cost, is not proven correct nor open-source. Hur et al. (2014) developed a different theory for slicing probabilistic programs. In their unpublished implementation, they observe significant runtime improvements if statements irrelevant to the return expression are removed from the program. In contrast, we slice programs with respect to the dependencies between sample statements in the program and there are no redundancies in the programs of our benchmark set. Castellan and Paquet (2019) present a theoretical foundation for incremental computation in the Metropolis Hastings algorithm by graphically modelling the data-dependencies between sample statements. However, they only consider a first-order probabilistic programs and did not implement or evaluate their approach.

Furthermore, there are many more optimisation opportunities to reduce the variance in variational inference like control variates (Paisley et al., 2012) or path-derivatives (Roeder et al., 2017). For a review on advances in variational inference see Zhang et al. (2018). Also, there are many methods to facilitate iterative implementations of SMC: WebPPL (Goodman and Stuhlmüller, 2014) and Anglican (Tolpin et al., 2015) use a continuation-passing-style transform; Turing (Ge et al., 2018) uses a tape of instructions; some PPLs like Birch (Murray and Schön, 2018) require the user to explicitly write the model in an iterative fashion.

Lastly, the probabilistic dependency structure of programs has also been exploited for applications different from optimising inference. Baah et al. (2008) leverage the structure for fault diagnosis and Bernstein et al. (2020) extract the factor graph of programs to automate model checking for models implemented in Stan. In general, static analysis for probabilistic programming is an active research field. For instance, Lee et al. (2019) developed a static analysis to verify so-called guide programs for stochastic variational inference in Pyro. Gorinova et al. (2021) developed an information flow type system and program transformation to automatically marginalise out discrete variables in Stan. Another example is the extension of the Stan compiler with a pedantic mode that statically catches common pitfalls in probabilistic programming (Bernstein, 2023). Recently, a framework for statically analysing probabilistic programs in a language-agnostic manner was presented (Böck et al., 2024). For a 2019 survey on static analysis for probabilistic programming see Bernstein (2019).

7. Conclusion

In this work, we addressed the open question to what extent a probabilistic program with user-labelled sample statements and while loops can be represented graphically. We presented a generalisation of existing formal semantics to support these language features. We developed a sound static analysis that approximates the dependency structure of random variables in the program and used this structure to factorise the implicitly defined program density. The factorisation is equivalent to the known Bayesian network factorisation in a restricted case, but novel for programs with dynamic addresses and loops. Lastly, we demonstrated that by generating a sub-program for each factor based on the static dependency structure, we can 1. speed-up a single-site Metropolis Hastings algorithm; 2. reduce gradient variance in variational inference; 3. facilitate an iterative implementation of sequential Monte Carlo. This work is to first to make these optimisations statically available for the considered program class. We proved correctness and empirically showed that our approach performs comparably or better than existing methods.

8. Data-Availability Statement

The source code to reproduce the experiments of Section 5 is made available at the time of review (Böck, 2025), will be submitted to Artifact Evaluation, and will be published after the review period. The shared repository includes the implementation of the static provenance analysis (Section 3.2) for a research PPL, the generation of sub-programs according to the dependency structure (Section 5), as well as a standard and optimised implementations of the LMH, BBVI, and SMC algorithms. It also includes the 16 benchmark programs (Section 5.2.1) written in our research PPL, with and without combinators in Gen, in recursive fashion in WebPPL, and in Pyro.

Appendix A Control-Flow Graph Construction

Definition 0.

A CFG is a tuple $(N_{\textnormal{start}},N_{\textnormal{end}},\mathcal{N},\mathcal{E})$ comprised of start and end node, a set of branch, join, and assign nodes, $\mathcal{N}$, and a set of edges, $\mathcal{E}\subseteq\mathcal{N}^{*}\times\mathcal{N}^{*}$, where $\mathcal{N}^{*}=\mathcal{N}\cup\{N_{\textnormal{start}},N_{\textnormal{end}}\}$.

• •

The CFG of a skip-statement, ${S=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{skip}}}}}}$, consists of only start and end node:

G=\big{(}N_{\textnormal{start}},N_{\textnormal{end}},\emptyset,\{(N_{\textnormal{start}},N_{\textnormal{end}})\}\big{)}

• •

The CFG of an assignment, $x=E$, or sample statement ${x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},f(E_{1},\dots,E_{n}))$, is a sequence of start, assign, and end node:

[TABLE]

where $A=\textnormal{Assign}(x=\dots)$ is the corresponding assign node.

• •

The CFG for a sequence of statements, $S=S_{1};S_{2}$, is recursively defined by ”stitching together” the CFGs $G_{1}=(N_{\textnormal{start}}^{1},N_{\textnormal{end}}^{1},\mathcal{N}^{1},\mathcal{E}^{1})$ and $G_{2}=(N_{\textnormal{start}}^{2},N_{\textnormal{end}}^{2},\mathcal{N}^{2},\mathcal{E}^{2})$ of $S_{1}$ and $S_{2}$:

[TABLE]

where $N^{1}_{e}$ is the predecessor node of $N_{\textnormal{end}}^{1}$ in $G_{1}$, $(N^{1}_{e},N_{\textnormal{end}}^{1})\in\mathcal{E}^{1}$, and $N^{2}_{s}$ is the successor node of $N_{\textnormal{start}}^{2}$ in $G_{2}$, $(N_{\textnormal{start}}^{2},N^{2}_{s})\in\mathcal{E}_{2}$.

• •

The CFG of an if statement, ${{{S=(\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{if}}}}}}\;E\;\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{then}}}}}}\;S_{1}\,\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{else}}}}}}\;S_{2})$, is also defined in terms of the CFGs of $S_{1}$ and $S_{2}$ together with a branch node $B=\textnormal{Branch}(E)$ and a join node $J$:

[TABLE]

The nodes $N^{1}_{s}$, $N^{1}_{e}$, $N^{2}_{s}$, $N^{2}_{e}$ are defined as in the last case. This construction assumes that $\mathcal{E}^{1}\neq\emptyset$ and $\mathcal{E}^{2}\neq\emptyset$. Otherwise, the branch would be a skip statement and we would connect $B$ with $J$ directly. The two predecessors of $J$ are denoted with $\textnormal{predcon}(J)$ and $\textnormal{predalt}(J)$ depending on the path.

• •

The CFG of ${{(\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{while}}}}}}\;E\;\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{do}}}}}}\;S)$ is given in terms of the sub-cfg $G_{S}$ of statement $S$.

[TABLE]

where $G_{S}=(N_{\textnormal{start}}^{S},N_{\textnormal{end}}^{S},\mathcal{N}^{S},\mathcal{E}^{S})$.

Again, if ${S=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{skip}}}}}}$, then we would have the edge $(B,B)$ instead of $(B,N^{S}_{s})$ and $(N^{S}_{e},B)$.

• •

The unrolled CFG of ${{(\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{while}}}}}}\;E\;\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{do}}}}}}\;S)$ is defined as follows:

[TABLE]

where $G_{S_{i}}=(N_{\textnormal{start}}^{S_{i}},N_{\textnormal{end}}^{S_{i}},\mathcal{N}^{S_{i}},\mathcal{E}^{S_{i}})$ is the $i$-th copy of sub-graph $G_{S}$. $N^{S_{i}}_{s}$ and $N^{S_{i}}_{e}$ are defined as in the cases before.

If ${S=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{skip}}}}}}$, then we would have the edges $(B_{i},B_{i+1})$ instead of the edges $(B_{i},N^{S_{i}}_{s})$ and $(N^{S_{i}}_{e},B_{i+1})$.

Appendix B Proofs of Sections 2-4

B.1. Proof of Proposition 2

Statement.* For all programs $S$ without while loops and corresponding control-flow graph $G$, it holds that for all traces $\mathtt{tr}$ and all initial states $\sigma_{0}$*

[TABLE]

Proof.

With structural induction we prove the more general equivalence:

[TABLE]

We write $\underset{S}{\overset{\mathtt{tr}}{\longrightarrow}}$ to indicate that the predecessor-successor relationship follow the CFG of program $S$.

For programs consisting of a single skip, assignment, or sample statement, the equality can be seen by directly comparing the operational semantic rules to the graph semantics.

For sequence $S=(S_{1};S_{2})$, by induction assumption we have

[TABLE]

By definition of the CFG of $S$, we see that

[TABLE]

In a similar way, we proof the result for if statements ${{{S=(\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{if}}}}}}\;E\;\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{then}}}}}}\;S_{1}\,\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{else}}}}}}\;S_{2})$ by considering the two cases $\sigma(E)=\mathtt{true}$ and $\sigma(E)=\mathtt{false}$.

[TABLE]

Depending on the case, we replace the dots with the path of $S_{i}$ without the start and end nodes. Note that in our notation in both cases $(\sigma,S_{i})\Downarrow^{\mathtt{tr}}\sigma^{\prime}$. ∎

B.2. Proof of Proposition 9

To prove soundness of Algorithm 1, we list properties of the provenance set that will become useful in subsequent proofs.

Lemma 0.

Let $N$ be a CFG node and $x$ a variable. prov has following properties:

• •

For all nodes $N^{\prime}$ it holds that

[TABLE]

• •

If $N^{\prime}\in\textnormal{RD}(N,x)$ is an assignment node, $N^{\prime}=\textnormal{Assign}(x=E)$, then

[TABLE]

• •

if $N^{\prime}\in\textnormal{RD}(N,x)$ is a sample node, ${N^{\prime}=\textnormal{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},\dots))$, then

[TABLE]

• •

For all $N^{\prime}\in\textnormal{RD}(N,x)$ it holds that for all branch parents $N_{\text{bp}}=\textnormal{Branch}(E_{N_{\text{bp}}})\in\textnormal{BP}(N^{\prime})$ we have $\textnormal{prov}(N_{\text{bp}},E_{N_{\text{bp}}})\subseteq\textnormal{prov}(N,x)$ and

[TABLE]

Proof.

If $(N^{\prime},y)$ is pushed into the queue, then $\textnormal{prov}(N^{\prime},y)\subseteq\textnormal{prov}(N,x)$ as the algorithm for computing $\textnormal{prov}(N^{\prime},y)$ starts with only $(N^{\prime},y)$ in the queue. With this fact the properties follow directly from the definition of the algorithm. ∎

Statement.* Let $G$ be the CFG of a program $S$ without while loop statements. For each node $N$ and variable $x$, there exists an evaluation function*

[TABLE]

such that for all traces $\mathtt{tr}$, initial program state $\sigma_{0}$ as in Definition 1, and execution sequence

[TABLE]

we have

[TABLE]

In the proof we also use the evaluation function of expressions $V^{N}_{E}\in[\mathcal{T}|_{\textnormal{prov}(N,E)}\to\mathcal{V}]$, where $V^{N}_{c}(\mathtt{tr})=c$ and $V^{N}_{g(E_{1},\dots,E_{n})}(\mathtt{tr})=g(V^{N}_{E_{1}}(\mathtt{tr}),\dots,V^{N}_{E_{n}}(\mathtt{tr})).$

Proof.

Step 1. Defining $V^{N}_{x}$ and proving $V^{N}_{x}\in[\mathcal{T}|_{\textnormal{prov}(N,x)}\to\mathcal{V}]$.

We begin by defining $V^{N}_{x}$ inductively by noting that every node $N\neq\text{START}$ has exactly one predecessor except for join nodes which have two predecessors. Further, the CFG $G$ is a finite directed acyclic graph with a single root node which makes mathematical induction possible.

Base case. For $((x_{i}\mapsto\mathtt{null},\textbf{p}\mapsto 1),\text{START})$, let $V^{\text{START}}_{\textbf{p}}(\mathtt{tr})=1$ and $V^{\text{START}}_{x}(\mathtt{tr})=\mathtt{null}$ for all variables $x$. We have $V^{\text{START}}_{x}\in[\mathcal{T}|_{\emptyset}\to\mathcal{V}]$.

Induction step.

Case 1. Let $N$ be a node with single predecessor $N^{\prime}$ ($N$ is not a join node). We define $V^{N}_{x}$ based on the node type of $N^{\prime}$, for which we make the induction assumption that the evaluation function $V^{N^{\prime}}_{x}\in[\mathcal{T}|_{\textnormal{prov}(N^{\prime},x)}\to\mathcal{V}]$ exists for all variables.

• •

Case 1.1. $N^{\prime}$ is an assignment node, $N^{\prime}=\textnormal{Assign}(z=E)$. Define

[TABLE]

By assumption $V^{N^{\prime}}_{y}\in[\mathcal{T}|_{\textnormal{prov}(N^{\prime},y)}\to\mathcal{V}]$ for all $y$ and thus, $V^{N^{\prime}}_{E}\in[\mathcal{T}|_{\textnormal{prov}(N^{\prime},E)}\to\mathcal{V}]$.

For $x\neq z$, $\textnormal{RD}(N,x)=\textnormal{RD}(N^{\prime},x)$ and by Eq. 11 we have $\textnormal{prov}(N,x)=\textnormal{prov}(N^{\prime},x)$, implying that $V^{N}_{x}\in[\mathcal{T}|_{\textnormal{prov}(N,x)}\to\mathcal{V}]$.

For $x=z$, $\textnormal{RD}(N,x)=\{N^{\prime}\}$ since $N^{\prime}$ is the single predecessor of $N$ and assigns $x$. By Eq. 12 $\textnormal{prov}(N^{\prime},E)\subseteq\textnormal{prov}(N,x)$, and by Eq. 4

[TABLE]

• •

Case 1.2. $N^{\prime}$ is a sample node, ${N^{\prime}=\textnormal{Assign}(z=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},f(E_{1},\dots,E_{n})))$. Define

[TABLE]

As before, for $x\neq z$, $V^{N}_{x}\in[\mathcal{T}|_{\textnormal{prov}(N,x)}\to\mathcal{V}]$.

For $x=z$, by definition $V^{N^{\prime}}_{E_{0}}(\mathtt{tr})\in\textnormal{addresses}(N^{\prime})$. Since $N^{\prime}\in\textnormal{RD}(N,x)=\{N^{\prime}\}$ it follows from Eq. 13 that $\textnormal{addresses}(N^{\prime})\subseteq\textnormal{prov}(N,x)$ and $\textnormal{prov}(N^{\prime},E_{0})\subseteq\textnormal{prov}(N,x)$. By Eq. 4

[TABLE]

• •

Case 1.3. $N^{\prime}$ is a branch or join node. Define $V^{N}_{y}=V^{N^{\prime}}_{y}$ for all variables $y$.

Case 2. Let $J$ be a join node with two predecessor nodes, $N^{\prime}_{1}=\text{predcons}(J)$, $N^{\prime}_{2}=\text{predalt}(J)$. Let $B=\text{Branch}(E)$ be the corresponding branching node, $(B,J)\in\textnormal{BranchJoin}(G)$.

For variable $x$, there are two cases:

1. (1)

There exists a reaching definition $N^{\prime}\in\textnormal{RD}(J,x)$ on a path $(B,\dots,N^{\prime},\dots,J)$. Define $V^{1}_{x}$ and $V^{2}_{x}$ depending on the node type of $N^{\prime}_{1}$ and $N^{\prime}_{2}$ as in step 1. The functions $V^{i}_{x}$ compute the value of $x$ after executing $N^{\prime}_{i}$. As before one can see that $V^{i}_{x}\in[\mathcal{T}|_{\textnormal{prov}(J,x)}\to\mathcal{V}]$. Let

[TABLE]

Since $B\in\textnormal{BP}(N^{\prime})$, by Eq. 14 $\textnormal{prov}(B,E)\subseteq\textnormal{prov}(J,x)$ and thus $V^{J}_{x}\in[\mathcal{T}|_{\textnormal{prov}(J,x)}\to\mathcal{V}]$. 2. (2)

All reaching definitions of $x$ (if there are any) are predecessors of $B$. Then, $\textnormal{RD}(J,x)=\textnormal{RD}(B,x)$ and by Eq. 11 $\textnormal{prov}(J,x)=\textnormal{prov}(B,x)$. Define $V^{J}_{x}\coloneqq V^{B}_{x}.$

Step 2: Proving that $\sigma_{i}(x)=V^{N_{i}}_{x}(\mathtt{tr})$.

Having defined the evaluation functions, we now have to prove that they indeed compute the correct values. For trace $\mathtt{tr}$ let the corresponding execution sequence be

[TABLE]

We will prove that $\sigma_{i}(x)=V^{N_{i}}_{x}(\mathtt{tr})$ by induction.

Base case. The base case immediately follows from the definition $\sigma_{0}(x)=V^{\text{START}}_{x}(\mathtt{tr})$.

Induction step. For transition $(\sigma_{i},N_{i})\overset{\mathtt{tr}}{\longrightarrow}(\sigma_{i+1},N_{i+1})$ the assumption is that $\sigma_{j}(x)=V^{N_{j}}_{x}(\mathtt{tr})$ holds for all variables $x$ and $j\leq i$.

If $N_{i+1}$ is not a join node, then $\sigma_{i+1}(x)=V^{N_{i+1}}_{x}(\mathtt{tr})$ follows directly from the CFG semantics and definition of $V^{N_{i+1}}_{x}$ in case 1 of step 1.

Lastly, we consider the case when $N_{i+1}$ is a join node with branching node $N_{j}=\text{Branch}(E)$, $(N_{j},N_{i+1})\in\text{BranchJoin}(G)$, for some $j\leq i$. Let $N^{\prime}_{1}=\text{predcons}(N_{i+1})$, $N^{\prime}_{2}=\text{predalt}(N_{i+1})$.

By the semantics of if statements, $V^{B}_{E}(\mathtt{tr})=\mathtt{true}$ iff $N_{i}=N^{\prime}_{1}$ and $V^{B}_{E}(\mathtt{tr})=\mathtt{false}$ iff $N_{i}=N^{\prime}_{2}$. If there is a node $N^{\prime}$ between $N_{j}$ and $N_{i+1}$ in the execution sequence that assigns $x$, $N^{\prime}\in\textnormal{RD}(N_{i+1},x)$, then by definitions of $V^{1}_{x}$, $V^{2}_{x}$, and $V^{N_{i+1}}_{x}$ we have

[TABLE]

If there is no such $N^{\prime}$, then

[TABLE]

∎

B.3. Proof of Proposition 1

Statement.* For all program $S$ with CFG $G$ and unrolled CFG $H$, it holds that for all traces $\mathtt{tr}$ and initial states $\sigma_{0}$*

[TABLE]

Proof.

We continue the proof of Proposition 2 for a while statement. Note the well-known equivalence for while statements in terms of the operational semantics:

[TABLE]

Thus, if the program terminates for $\mathtt{tr}$, it is equivalent to $S;\dots;S$ ($S$ repeated $n$ times), where $\sigma_{0}=\sigma$, $(\sigma_{i},S)\Downarrow^{\mathtt{tr}}\sigma_{i+1}$, $\sigma_{n}=\sigma^{\prime}$, such that $\sigma_{i}(E)=\mathtt{true}$ for $i<n$, $\sigma_{n}(E)=\mathtt{false}$.

In this case, it can be shown that the (unrolled) CFG of $S;\dots;S$ is also equivalent to the (unrolled) CFG of the while statement.

If the program does not terminate, then ${\bar{p}_{S}(\mathtt{tr},\sigma_{0})=\bar{p}_{G}(\mathtt{tr},\sigma_{0})=\bar{p}_{H}(\mathtt{tr},\sigma_{0})}=\textnormal{undefined}$. ∎

B.4. Proof of Lemma 3

Statement.* Let $G$ be the CFG and $H$ the unrolled CFG of program $S$. For each node $M\in H$ following equation holds.*

[TABLE]

Proof.

We will prove following Lemma by structural induction:

Lemma. For every path $\vv{M}=(M^{\prime},\dots,M)$ in the unrolled CFG $H$, there is a path in the CFG $G$, $\vv{N}=(N^{\prime},\dots,N)$, such that

[TABLE]

[TABLE]

*for all variables $x$, where $\textnormal{AssignNodes}(\vv{N},x)$ is the set of all CFG nodes that assign $x$ in path $\vv{N}$ and $\textnormal{AssignNodes}_{\textnormal{cfg}}(\vv{M},x)$ is the set of all assign nodes in path $\vv{M}$ mapped to $G$ with $\iota_{\textnormal{cfg}}$. *

With this lemma we can prove the statement:

Let $M^{\prime}\in\textnormal{RD}(M,x)$. Thus, there exists a path $\vv{M}=(M^{\prime},\dots,M)$, such that $M^{\prime}$ is the only assign node for $x$ before $M$. By the lemma, there is a path in the CFG $\vv{N}=(N^{\prime},\dots,N)$, such that $\iota_{\textnormal{cfg}}(M^{\prime})=N^{\prime}$ and $\iota_{\textnormal{cfg}}(M)=N$.

Further, since $\textnormal{AssignNodes}(\vv{N},x)=\textnormal{AssignNodes}_{\textnormal{cfg}}(\vv{M},x)\subseteq\{\iota_{\textnormal{cfg}}(M^{\prime}),\iota_{\textnormal{cfg}}(M)\}=\{N^{\prime},N\}$, the node $N^{\prime}$ is also the only node in the path $\vv{N}$ before $N$ that assigns $x$.

Thus, $\iota_{\textnormal{cfg}}(M^{\prime})=N^{\prime}\in\textnormal{RD}(N,x)=\textnormal{RD}(\iota_{\textnormal{cfg}}(M),x).$

Proof of Lemma. By structural induction.

For CFGs consisting only of a single assign node, the unrolled CFGs are identical and the statement follows.

For a sequence of two statements $(S_{1};S_{2})$, let $G$ be the CFG with sub-graphs $G_{1}$ and $G_{2}$ for statements $S_{1}$ and $S_{2}$ respectively. Let $H$ the unrolled CFG with sub-graphs $H_{1}$ and $H_{2}$. We consider the only interesting case where $M^{\prime}\in H_{1}$ and $M\in H_{2}$:

[TABLE]

For sub-paths $\vv{M_{1}}\coloneqq(M,\dots,M_{i},\textnormal{END})$ and $\vv{M_{2}}\coloneqq(\textnormal{START},M_{i+1},\dots,M)$ we apply the induction assumption to get two paths $\vv{N_{1}}=(N^{\prime},\dots,N_{1},\text{END})$, $\vv{N_{2}}=(\text{START},N_{2},\dots,N)$, with $\iota_{\textnormal{cfg}}(M^{\prime})=N^{\prime}$ and $\iota_{\textnormal{cfg}}(M)=N$. The corresponding CFG path in $G$ is $\vv{N}\coloneqq(N^{\prime},\dots,N_{1},N_{2},\dots,N)$ as

[TABLE]

For an if statement, let $G$ be the CFG with branch node $\tilde{B}$, join node $\tilde{J}$, and sub-graphs $G_{1}$, $G_{2}$. Let $H$ be the unrolled CFG with corresponding $B$, $J$, $H_{1}$, and $H_{2}$. We only show the statement for the case

[TABLE]

Again, we apply the induction assumption to $\vv{M_{i}}\coloneqq(\textnormal{START},M^{\prime},\dots,M,\textnormal{END})$ and get a path in $G_{i}$: $\vv{N_{i}}=(\textnormal{START},N^{\prime},\dots,N,\textnormal{END})$. Since $\iota_{\textnormal{cfg}}(B)=\tilde{B}$, $\iota_{\textnormal{cfg}}(J)=\tilde{J}$, which are not assign nodes, and $\textnormal{AssignNodes}(\vv{N_{i}})=\textnormal{AssignNodes}_{\textnormal{cfg}}(\vv{M_{i}})$, the statement follows for the path $(\tilde{B},N^{\prime},\dots,N,\tilde{J})$.

Lastly, for a while loop, let $G$ be the CFG with branch node $\tilde{B}$, join node $\tilde{J}$, and sub-graph $G_{S}$. Let $H$ be the unrolled CFG with $B_{i}$ branch nodes, $J$ join node, and $H_{i}$ sub-graphs. For the sake of brevity we consider three cases:

• •

For $M^{\prime}\in H_{i}$ and $M\in H_{j}$

[TABLE]

we apply the induction assumption to the sub-paths and construct the path

[TABLE]

where $\tilde{B}$ and $B_{i}$ are do not contribute to $\textnormal{AssignNodes}(\vv{N})$ or $\textnormal{AssignNodes}(\vv{M})$.

• •

For $M^{\prime}\in H_{i}$ and $J$

[TABLE]

we construct the path

[TABLE]

• •

For the path $\vv{M}=(B_{l},J)$ the path in the CFG is $\vv{N}=(\tilde{B},\tilde{J})$.

∎

Statement.* Let $G$ be the CFG and $H$ the unrolled CFG of program $S$. For each node $M\in H$ the following equation holds.*

[TABLE]

Proof.

We prove this statement by structural induction.

For CFGs only consisting of a single assign node, the unrolled CFGs are identical, there are no branch parents, and the statement follows.

The sequencing of two statements $(S_{1};S_{2})$ does not change the branch parents of any node.

For an if statement, let $G$ be the CFG with branch node $\tilde{B}$, join node $\tilde{J}$, and sub-graphs $G_{1}$, $G_{2}$. Let $H$ be the unrolled CFG with corresponding $B$, $J$, $H_{1}$, and $H_{2}$. Every node $M$ in $H_{i}$ has branch parents $\left(\textnormal{BP}(M)\cap H_{i}\right)\cup\{B\}$.

By induction assumption $\textnormal{BP}_{\textnormal{cfg}}(M)\cap G_{i}=\textnormal{BP}(\iota_{\textnormal{cfg}}(M))\cap G_{i}$ and thus

[TABLE]

Nodes $B$, $\tilde{B}$, $J$, and $\tilde{J}$ do not have branch parents.

Lastly, for a while loop, let $G$ be the CFG with branch node $\tilde{B}$, join node $\tilde{J}$, and sub-graph $G_{S}$. Let $H$ be the unrolled CFG with $B_{i}$ branch nodes, $J$ join node, and $H_{i}$ sub-graphs.

Every node $M$ in $H_{i}$ has branch parents $\left(\textnormal{BP}(M)\cap H_{i}\right)\cup\{B_{1},\dots,B_{i}\}$.

The CFG node $\iota_{\textnormal{cfg}}(M)\in G_{S}$ has branch parents $\left(\textnormal{BP}(\iota_{\textnormal{cfg}}(M))\cap G_{S}\right)\cup\{\tilde{B}\}$.

Since $\iota_{\textnormal{cfg}}(B_{j})=\tilde{B}$ and by induction assumption $\textnormal{BP}_{\textnormal{cfg}}(M)\cap G_{S}=\textnormal{BP}(\iota_{\textnormal{cfg}}(M))\cap G_{S}$, we get

[TABLE]

∎

Statement.* Let $G$ be the CFG and $H$ the unrolled CFG of program $S$. For each node $M\in H$ the following equation holds.*

[TABLE]

Proof.

Follows directly from the two equations $\{\iota_{\textnormal{cfg}}(M^{\prime})\colon M^{\prime}\in\textnormal{RD}(M,x)\}\subseteq\textnormal{RD}(\iota_{\textnormal{cfg}}(M),x)$ and $\{\iota_{\textnormal{cfg}}(M^{\prime})\colon M^{\prime}\in\textnormal{BP}(M)\}=\textnormal{BP}(\iota_{\textnormal{cfg}}(M))$, since $\textnormal{prov}(M,x)$ is defined in terms of RD and BP. ∎

B.5. Proof of Proposition 4

Statement.* Let $H$ be the unrolled CFG of a program $S$. For each node $M$ and variable $x$, there exists an evaluation function*

[TABLE]

such that for all traces $\mathtt{tr}$, initial program state $\sigma_{0}$ as in Definition 1, and execution sequence

[TABLE]

we have

[TABLE]

Proof.

First, we have to justify the use of mathematical induction to prove properties on a potentially infinite graph. For $M_{1},M_{2}\in H$ define the relation $M_{1}\prec M_{2}\Leftrightarrow M_{1}\text{ is parent of }M_{2}$ (i.e. there is an edge from $M_{1}$ to $M_{2}$ in $H$). By construction of $H$, there are no infinite descending chains, $M_{i}\in H$ for $i\in\mathbb{N}$ such that $M_{i+1}\prec M_{i}$. This makes $\prec$ a well-founded relation, i.e. every non-empty subset has a minimal element. As a result, we may apply induction with respect to this relation (as has been done in the proof of Proposition 9) for graphs with infinite nodes. This proof principle is called well-founded or Noetherian induction.

To complete the proof, we only need to consider one additional case:

Case 3. Let $J$ be a join node of a while loop with the parents $B_{i}=\textnormal{Branch}(E)$, $i\in\mathbb{N}$. The function $V^{B_{i}}_{x}\in[\mathcal{T}|_{\textnormal{prov}(B_{i},x)}\to\mathcal{V}]$ exists by induction assumption and computes the value of $x$ before executing the while loop body for the $i$-th time. If there is a reaching definition $M^{\prime}\in\textnormal{RD}(J,x)$ on a path $(B_{i},\dots,M^{\prime},\dots,J)$, define

[TABLE]

Further, as $M^{\prime}$ belongs to the sub-CFG of the while loop body, we have that for every $B_{i}$ there is a $M_{i}^{\prime}\in\textnormal{RD}(J,x)$ on the path from $B_{i}$ to $J$. Since $B_{i}\in\textnormal{BP}(M_{i}^{\prime})$, by Equation 14 $\textnormal{prov}(B_{i},E)\subseteq\textnormal{prov}(J,x)$ for all $i$ and thus $V^{J}_{x}\in[\mathcal{T}|_{\textnormal{prov}(J,x)}\to\mathcal{V}]$. If there is no reaching definition in the while loop body sub-CFG, then $\textnormal{prov}(J,x)=\textnormal{prov}(B_{1},x)$ and define $V^{J}_{x}=V^{B_{1}}_{x}$.

∎

B.6. Proof of Theorem 5

Statement.* Let $G$ be the CFG for a program $S$. Let $N_{1},\dots,N_{K}$ be all sample nodes in $G$. For each sample node ${N_{k}=\textnormal{Assign}(x_{k}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{k},f^{k}(E_{1}^{k},\dots,E_{n_{k}}^{k})))$, let*

[TABLE]

Then, there exist functions $p_{k}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}]$ such that for all $\mathtt{tr}\in\mathcal{T}$, if $p_{S}(\mathtt{tr})\neq\textnormal{undefined}$, then

[TABLE]

Proof.

Let $H$ be the unrolled CFG of program $S$. Denote the (countably many) sample nodes in $H$ with ${M_{j}=\textnormal{Assign}(x^{j}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{j},f^{j}(E_{1}^{j},\dots,E_{n_{j}}^{j})))$. Like in the proof of Theorem 10, define for each $M_{j}$, $b_{j}(\mathtt{tr})\coloneqq\bigwedge_{M^{\prime}\in\textnormal{BP}(M_{j})}t^{j}_{M^{\prime}}(V^{M^{\prime}}_{{\textnormal{condexp}(M^{\prime})}}(\mathtt{tr}))$ where

[TABLE]

Again, $b_{j}(\mathtt{tr})=\mathtt{true}$ if $M_{j}$ is in the execution sequence for $\mathtt{tr}$ else $\mathtt{false}$. Define

[TABLE]

We have $\tilde{p}_{j}\in[\mathcal{T}|_{\tilde{A}_{j}}\to\mathbb{R}_{\geq 0}]$, for

[TABLE]

Finally, we group the nodes $M_{j}$ by their corresponding CFG node.

[TABLE]

If $\iota_{\textnormal{cfg}}(M_{j})=N_{k}$, then $\tilde{A}_{j}\subseteq A_{k}$, since for $i=0,\dots,n$ we have $E_{i}^{k}=E_{i}^{j}$ and by Eq. 8

[TABLE]

and by Eq. 7 it holds that

[TABLE]

such that

[TABLE]

Thus, for all $j$, $\iota_{\textnormal{cfg}}(M_{j})=N_{k}$, we have $\tilde{p}_{j}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}]$ which implies $p_{k}\in[\mathcal{T}|_{A_{k}}\to\mathbb{R}_{\geq 0}].$ Again, we know that for the end-state $\sigma_{l}(\textbf{p})=V^{\text{END}}_{\textbf{p}}(\mathtt{tr})=p_{H}(\mathtt{tr})$ and since $b_{j}(\mathtt{tr})=\mathtt{true}$ if and only if $M_{j}$ is in the execution sequence for $\mathtt{tr}$, we have $p_{H}(\mathtt{tr})=\prod_{k=1}^{K}p_{k}(\mathtt{tr})$. ∎

Appendix C Towards a Measure-theoretic Interpretation of the Operational Semantics

In this section, we will briefly discuss how one may construct a measure space on traces such that the function $\mathtt{tr}\mapsto p_{S}(\mathtt{tr})$ is indeed a Radon-Nikodym derivative of a measure on this space. This construction is related to denotational semantics based on trace types (Lew et al., 2019, 2023). We first define a measure space for each value type:

[TABLE]

These measure spaces model the support of common distributions, but could also be extended by building sum and product spaces. They are comprised of the typical Dirac measure, power sets, counting measure, Borel sets, and Lebesgue measures.

Each trace $\mathtt{tr}$ can be assigned an unique type function $\tau_{\mathtt{tr}}\colon\mathbf{Strings}\to\mathbf{Types}$ such that for all $\alpha\in\mathbf{Strings}$ we have $\mathtt{tr}(\alpha)\in\tau_{\mathtt{tr}}(\alpha)$. The types are assumed to be disjoint (e.g. $\mathbb{Z}\cap\mathbb{R}=\emptyset$). Let

[TABLE]

be the set of all type functions that correspond to traces with a finite number of non-$\mathtt{null}$ values. We enumerate the addresses, $\mathbf{Strings}=\{\alpha_{i}\colon i\in\mathbb{N}\}$, and define a measureable space for each type function, $s\colon\mathbf{Strings}\to\mathbf{Types}$, $\mathscr{M}_{s}=(M_{s},\Sigma_{s})=\bigotimes_{i\in\mathbb{N}\colon s(\alpha_{i})\neq\{\mathtt{null}\}}\mathcal{M}_{s(\alpha_{i})}$. If $s\in\mathcal{S}_{<\infty}$, then $\mathscr{M}_{s}$ is a finite dimensional product space with product measure $\nu_{s}$. For the set of traces $\mathcal{T}$, we define the $\sigma$-algebra $\Sigma_{\mathcal{T}}$ with $A\in\Sigma_{\mathcal{T}}\Leftrightarrow\forall s\colon\mathbf{Strings}\to\mathbf{Types}\colon\pi_{s}(\{\mathtt{tr}\in A\colon\tau_{\mathtt{tr}}=s\})\in\Sigma_{s},$ where $\pi_{s}\colon\mathcal{T}\to M_{s}$ is the canonical projection. The reference measure on $(\mathcal{T},\Sigma_{\mathcal{T}})$ is

[TABLE]

which effectively is a measure of traces with a finite number of non-$\mathtt{null}$ values.

If a primitive distribution $f$ with $\text{support}(f)\subseteq\tau_{f}$ and arguments $v_{1}\dots,v_{n}$ corresponds to the measure $\mu_{f,v_{1},\dots,v_{n}}$ with Radon-Nikodym derivative ${\mathtt{d}\mu_{f,v_{1},\dots,v_{n}}}/{\mathtt{d}\nu_{\tau_{f}}}$, we define

[TABLE]

The measure-theoretic interpretation of a program is the measure $\mu_{S}(A)=\int_{\mathcal{T}_{\textnormal{valid}}\cap A}p_{S}\,\mathtt{d}\nu$, where

[TABLE]

is the set of all traces with well-defined density such that changing a value at any address to $\mathtt{null}$ leads to an undefined density. That is, for $\mathtt{tr}\in\mathcal{T}_{\textnormal{valid}}$ all addresses that are not used during the execution of the program are set to $\mathtt{null}$ and all address that are used are set to non-$\mathtt{null}$ values. Such a trace necessarily has a finite number of non-$\mathtt{null}$ values as the density for non-termination is undefined which is the only possibility for executing an infinite number of sample statements. This condition is important, because the countably infinite sum in definition of $\nu$ has better convergence properties on $\mathcal{T}_{\textnormal{valid}}$. Finally, $p_{S}$ the Radon-Nikodym derivative of $\mu_{S}$ with respect to $\nu$ on $\mathcal{T}_{\textnormal{valid}}$.

We emphasize that this interpretation demands a careful analysis of the measurability of $\mathcal{T}_{\textnormal{valid}}$ and $p_{S}$ which is beyond the scope of this work. This requires additional assumptions about the program, at the very least that the primitives $g$ are measurable. Such measurability analysis may be done along the lines of Borgström et al. (2016) and Lee et al. (2019).

Note that this construction can handle programs where random variables are of mixed-type like the variable $\mathtt{y}$ in the example below:

⬇

x = sample("x", Bernoulli(p))

if x == 1 then

y = sample("y", Bernoulli(0.5))

else

y = sample("y", Normal(0.,1.))

We list the density according to the semantics for different example trace of different type:

[TABLE]

C.1. Remark on Conditioning

We have defined density semantics such that $p_{S}$ denotes the joint probability distribution over all variables in the model in the sense described above. In our construction, we do not differentiate between latent and observed variables. To incorporate observed data, we can define the trace $\mathtt{obs\_tr}$, where addresses are mapped to the observed data points. All addresses that do not correspond to observed data map to $\mathtt{null}$. Loosely, the function $\mathtt{tr}\mapsto p_{S}(\mathtt{tr}\oplus\mathtt{obs\_tr})$ then corresponds to the unnormalised posterior density, where we merge the trace $\mathtt{tr}$, encapsulating latent variables, with the observed trace $\mathtt{obs\_tr}$ like below

[TABLE]

This is only a loose correspondence, because it is not clear whether this function can be normalized by integrating over all observed variables. It is beyond the scope of this work to precisely define conditioning semantics. For work discussing similar issues, see Lee et al. (2019).

Appendix D Proofs of Section 5

D.1. Sliced CFG Construction

First, we formally define the CFG slicing with respect to the factorisation of Theorem 5. Here, we use provnode instead of prov to make the dependencies between CFG nodes explicit (see Section 3.2).

Definition 0.

*Let $G=(N_{\textnormal{start}},N_{\textnormal{end}},\mathcal{N},\mathcal{E})$ be a CFG with $K$ sample nodes, which we enumerate with ${N_{j}=\textnormal{Assign}(x_{j}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{j},f^{j}(E_{1}^{j},\dots,E_{n_{j}}^{j})))$ for $j=1,\dots,K$.

For each $k\in\{1,\dots,K\}$, we define the set of indices of sample nodes that the factor $p_{k}$ depends on:*

[TABLE]

Next, we define the set of nodes that lie on path from $N_{k}$ to some of its dependent node $N_{j}$:

[TABLE]

We transform all sample nodes that are in $\mathcal{N}_{k}^{\prime}$, but do not depend on $N_{k}$ to a read statement, which in contrast to sample statements do not contribute to the density: $\mathcal{N}_{k}\coloneqq\{\psi(N)\colon N\in\mathcal{N}_{k}^{\prime}\}$, where

[TABLE]

The semantic of read statements is given below. Note that we have to keep all $E_{i}$ in the statement to check their values against $\mathtt{null}$ like in Equation 1. We omit this technicality in the main text.

[TABLE]

We restrict the set of edges to the subset of nodes $\mathcal{N}_{k}$:

[TABLE]

We connect the new start node $N^{k}_{\textnormal{start}}$ to $N_{k}$ and the new end node $N^{k}_{\textnormal{end}}$ to all $N_{j}$, $j\in\mathcal{J}_{k}$, that do not have a successor in $\mathcal{E}_{k}^{\prime}$:

[TABLE]

Finally, define the CFG sliced with respect to $N_{k}$ as

[TABLE]

Note that in this construction, we allow that the end node $N^{k}_{\textnormal{end}}$ has multiple predecessors. This does not change the CFG semantics. Next, we rewrite Lemma 3 and Lemma 1 in terms of provnode and prove two technical lemmas needed to establish the correctness result in the next section.

Lemma 0.

Let $M$ be a CFG node and $x$ a variable. provnode has following properties:

• •

If $M^{\prime}\in\textnormal{RD}(M,x)$ is an assignment node, $M^{\prime}=\textnormal{Assign}(x=E)$, then

[TABLE]

• •

if $M^{\prime}\in\textnormal{RD}(M,x)$ is a sample node, ${M^{\prime}=\textnormal{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},\dots))$, then

[TABLE]

• •

For all $M^{\prime}\in\textnormal{RD}(M,x)$ it holds that for all branch parents $N_{\text{bp}}=\textnormal{Branch}(E_{M_{\text{bp}}})\in\textnormal{BP}(M^{\prime})$

[TABLE]

• •

We have following relationship between provnode on CFG $G$ and on unrolled CFG $H$:

[TABLE]

Lemma 0.

Let $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ be traces that only differ at address $\alpha$, i.e. $\mathtt{tr}_{1}(\alpha)\neq\mathtt{tr}_{2}(\alpha)$ and $\forall\beta\neq\alpha\colon\mathtt{tr}_{1}(\beta)=\mathtt{tr}_{2}(\beta).$ Further, let $M_{\alpha}$ bet the only sample node in the execution sequences of $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ whose address expression evaluates to $\alpha$. If CFG node $M$ appears in both execution sequences and variable $x$,

[TABLE]

*then the inequality $\sigma_{1}(x)=V^{M}_{x}(\mathtt{tr}_{1})\neq V^{M}_{x}(\mathtt{tr}_{2})=\sigma_{2}(x)$ implies $M_{\alpha}\in\textnormal{provnode}(M,x)$.

Or equivalently, $M_{\alpha}\notin\textnormal{provnode}(M,x)$ implies $\sigma_{1}(x)=V^{M}_{x}(\mathtt{tr}_{1})=V^{M}_{x}(\mathtt{tr}_{2})=\sigma_{2}(x)$.*

Proof.

By induction on the position in the execution sequence. The claim is trivially true for the start node. Let $M$ be a node in both execution sequences with $V^{M}_{x}(\mathtt{tr}_{1})\neq V^{M}_{x}(\mathtt{tr}_{2})$. The induction assumption states that for all nodes $M^{\prime}$ that appear sooner than $M$ in any of the execution sequences, the claim holds.

Let $M^{\prime}_{1}\in\textnormal{RD}(M,x)\neq\emptyset$ be the assign node that last writes $x$ in the execution sequence of $\mathtt{tr}_{1}$ before evaluating $M$. Similarly, let $M^{\prime}_{2}\in\textnormal{RD}(M,x)\neq\emptyset$ be the corresponding assign node in the execution sequences of $\mathtt{tr}_{2}$.

• •

Case 1. $M^{\prime}_{1}\neq M^{\prime}_{2}$: By definition of reaching definitions, both execution sequences must contain a different path from start node to $M$. This can only happen if there exists one branch parent $B\in\textnormal{BP}(M^{\prime}_{1})\cap\textnormal{BP}(M^{\prime}_{2})$ in both execution sequences that evaluates differently, $V^{B}_{E}(\mathtt{tr}_{1})\neq V^{B}_{E}(\mathtt{tr}_{2})$ for $E=\textnormal{condexp}(B)$. By induction assumption, $M_{\alpha}\in\textnormal{provnode}(B,E)$ and Equation 21 implies $M_{\alpha}\in\textnormal{provnode}(M,x)$.

• •

Case 2. $M^{\prime}\coloneqq M^{\prime}_{1}=M^{\prime}_{2}=\textnormal{Assign}(x=E)$: It holds that $V^{M^{\prime}}_{E}(\mathtt{tr}_{1})=V^{M}_{x}(\mathtt{tr}_{1})\neq V^{M}_{x}(\mathtt{tr}_{2})=V^{M^{\prime}}_{E}(\mathtt{tr}_{2})$. By induction assumption, we have $M_{\alpha}\in\textnormal{provnode}(M^{\prime},E)$ and Equation 19 implies $M_{\alpha}\in\textnormal{provnode}(M,x)$.

• •

Case 3. ${M^{\prime}\coloneqq M^{\prime}_{1}=M^{\prime}_{2}=\textnormal{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},\dots))$: It holds that $\mathtt{tr}_{1}(V^{M^{\prime}}_{E_{0}}(\mathtt{tr}_{1}))=V^{M}_{x}(\mathtt{tr}_{1})\neq V^{M}_{x}(\mathtt{tr}_{2})=\mathtt{tr}_{2}(V^{M^{\prime}}_{E_{0}}(\mathtt{tr}_{2}))$. If $V^{M^{\prime}}_{E_{0}}(\mathtt{tr}_{1})=\alpha$ or $V^{M^{\prime}}_{E_{0}}(\mathtt{tr}_{2})=\alpha$, then $M^{\prime}=M_{\alpha}$, since $M_{\alpha}$ is the only sample node whose address expression evaluates to $\alpha$ and we have $M_{\alpha}\in\textnormal{provnode}(M,x)$ by Equation 20. Otherwise, since the traces only differ at $\alpha$, the inequality can only hold if $V^{M^{\prime}}_{E_{0}}(\mathtt{tr}_{1})\neq V^{M^{\prime}}_{E_{0}}(\mathtt{tr}_{2})$. But then as in Case 2, by induction assumption $M_{\alpha}\in\textnormal{provnode}(M^{\prime},E_{0})$ and by Equation 20 we have $M_{\alpha}\in\textnormal{provnode}(M,x)$.

Lastly, if there are no such $M^{\prime}_{1}$ and $M^{\prime}_{2}$ in both sequences, then $V^{M}_{x}(\mathtt{tr}_{1})=\mathtt{null}=V^{M}_{x}(\mathtt{tr}_{2})$, a contradiction. If there is only $M^{\prime}_{1}$ or only $M^{\prime}_{2}$, then the claim follows as in Case 1. ∎

Lemma 0.

Let $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ be traces that only differ at address $\alpha$, i.e. $\mathtt{tr}_{1}(\alpha)\neq\mathtt{tr}_{2}(\alpha)$ and $\forall\beta\neq\alpha\colon\mathtt{tr}_{1}(\beta)=\mathtt{tr}_{2}(\beta).$ Further, let $M_{\alpha}$ be the only sample node in the execution sequences whose address expression evaluates to $\alpha$. If CFG node $M$ is in one but not both of the execution sequences of $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$, then there exists a branch parent $B\in\textnormal{BP}(M)$ such that $M_{\alpha}\in\textnormal{provnode}(B,\textnormal{condexp}(B))$.

Proof.

As in the proof of Theorem 5, $b(\mathtt{tr})=\bigwedge_{B\in\textnormal{BP}(M)}t_{B}(V^{B}_{\textnormal{condexp}(B)}(\mathtt{tr}))$ is $\mathtt{true}$ if $M$ is in the execution sequence of $\mathtt{tr}$ else $\mathtt{false}$. Since $M$ is not in both execution sequences we have $b(\mathtt{tr}_{1})\neq b(\mathtt{tr}_{2})$. Thus, one of the branch parents has to evaluate differently, i.e. there is a branch parent $B\in\textnormal{BP}(M)$ in both execution sequences such that $V^{B}_{\textnormal{condexp}(B)}(\mathtt{tr}_{1})\neq V^{B}_{\textnormal{condexp}(B)}(\mathtt{tr}_{2})$.

Since, $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ differ only at $\alpha$ and $M_{\alpha}$ is the only sample node in the execution sequences whose address expression evaluates to $\alpha$, by Lemma 3, $M_{\alpha}\in\textnormal{provnode}(B,\textnormal{condexp}(B))$. ∎

D.2. Proof of Theorem 1 - Correctness of Slicing

Now, we are ready to prove the main correctness result of the sliced CFG construction.

Statement.* Let $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ be traces that only differ at address $\alpha$, i.e. $\mathtt{tr}_{1}(\alpha)\neq\mathtt{tr}_{2}(\alpha)$ and $\forall\beta\neq\alpha\colon\mathtt{tr}_{1}(\beta)=\mathtt{tr}_{2}(\beta).$ Let the execution sequences in the unrolled CFG $H$ be equal up to sample node ${M_{\alpha}=\textnormal{Assign}(x=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0},f(E_{1},\dots,E_{n})))$ and state $\sigma$, such that $\alpha=\sigma(E_{0})$:*

[TABLE]

If $M_{\alpha}$ is the only sample node in the execution sequences whose address expression evaluates to $\alpha$, then for the sliced CFG $G_{k}$ of sample node $N_{k}=\iota_{\textnormal{cfg}}(M_{\alpha})$, it holds that

[TABLE]

Proof.

Let $M_{1}^{\prime}$ and $M_{2}^{\prime}$ be the last sample nodes in the sequences, such that $\iota_{\textnormal{cfg}}(M_{1}^{\prime})\in G_{k}$ and $\iota_{\textnormal{cfg}}(M_{2}^{\prime})\in G_{k}$. By definition of $G_{k}$, the execution sequences for $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ in the unrolled version $H_{k}$ of $G_{k}$ become

[TABLE]

[TABLE]

for some states $\sigma^{\prime}_{1}$, $\sigma^{\prime}_{2},\sigma^{\prime\prime}_{1}$, $\sigma^{\prime\prime}_{2}$. Note that technically we also have $(\sigma^{\prime}_{i},M^{\prime}_{i})\underset{H_{k}}{\overset{\mathtt{tr}_{i}}{\longrightarrow}}(\sigma_{i}^{\prime},\text{END})$.

To prove the statement, we have to show that the contribution to the density of sample nodes $M$ that come after $M_{1}^{\prime}$ and $M_{2}^{\prime}$ in the original execution sequences (23) or (24) and sample nodes $M$ that come between $M_{\alpha}$ and $M_{1}^{\prime}$ or $M_{2}^{\prime}$ but were transformed to read statements in $G_{k}$, see (18), is equal.

By Theorem 5, the contribution to the density of sample node $M$ with $N_{j}=\iota_{\textnormal{cfg}}(M)$, i.e. ${M=\textnormal{Assign}(x_{j}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{j},f^{j}(E_{1}^{j},\dots,E_{n_{j}}^{j})))$, has provenance set

[TABLE]

or equivalently expressed in terms of sample nodes

[TABLE]

such that by Lemma 2

[TABLE]

If $M$ comes after $M^{\prime}_{1}$ or $M^{\prime}_{2}$ in (23) or (24), then $\iota_{\textnormal{cfg}}(M)=N_{j}\notin G_{k}$, and by definition of $G_{k}$, we have $\iota_{\textnormal{cfg}}(M_{\alpha})=N_{k}\notin\mathcal{A}_{j}$ and thus also $M_{\alpha}\notin\tilde{\mathcal{A}}_{j}$. If $M$ comes between $M_{\alpha}$ and $M^{\prime}_{1}$ or $M^{\prime}_{2}$, but was transformed to a read statement in $G_{k}$, then by definition of $G_{k}$, also $N_{k}\notin\mathcal{A}_{j}$ and $M_{\alpha}\notin\tilde{\mathcal{A}}_{j}$.

In both cases, by Lemma 4, $M$ has to be in both original execution sequences, since otherwise we would have the contradiction $M_{\alpha}\in\textstyle\bigcup_{B\in\textnormal{BP}(M)}\textnormal{provnode}(B,\textnormal{condexp}(B))$. Let $(\sigma_{1},M)$ denote the appearance in (23) and $(\sigma_{2},M)$ the appearance in (24). Since $M_{\alpha}\notin\textstyle\bigcup_{i=0}^{n_{j}}\textnormal{provnode}(M,E_{i}^{j})$, by Lemma 3, we have $\sigma_{1}(y)=\sigma_{2}(y)$ for all variables $y\in\textstyle\bigcup_{i=0}^{n_{j}}\textnormal{vars}(E_{i}^{j})$. Since it also holds that $\alpha\neq\sigma_{1}(E_{0})=\sigma_{2}(E_{0})$, we have $\mathtt{tr}_{1}(\sigma_{1}(E_{0}))=\mathtt{tr}_{2}(\sigma_{2}(E_{0}))$, which means that $M$ contributes the same value to the density in both execution sequences, compare Equation 2. ∎

D.3. Proof of Corollary 2 - Correctness of Slicing (cont.)

Statement.* Assume that for each trace $\mathtt{tr}$ there is at most one sample node $M_{\alpha}$ in its execution sequences such that the address expression of $M_{\alpha}$ evaluates to $\alpha$. For initial program state $\sigma_{0}$ as in Definition 1, define*

[TABLE]

*Then, $\Delta_{\alpha}(\mathtt{tr})\coloneqq\log p_{G}(\mathtt{tr})-\log p_{\alpha}(\mathtt{tr})$ is independent of address $\alpha$, $\Delta_{\alpha}\in[\mathcal{T}|_{\mathbf{Strings}\setminus\{\alpha\}}\to\mathbb{R}_{\geq 0}].$ *

Proof.

By definition $\Delta_{\alpha}\in[\mathcal{T}|_{\mathbf{Strings}\setminus\{\alpha\}}\to\mathbb{R}_{\geq 0}]$, if for two traces $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ it holds that

[TABLE]

Let $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ be traces such that $\mathtt{tr}_{1}(\alpha)\neq\mathtt{tr}_{2}(\alpha)$ and $\forall\beta\neq\alpha\colon\mathtt{tr}_{1}(\beta)=\mathtt{tr}_{2}(\beta).$ The execution sequences of $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ are equal up to $(\sigma,M_{\alpha})$, where $M_{\alpha}$ is the first sample node whose address expression evaluates to $\alpha$. If there is no such sample node, then the execution sequences are completely identical. In this case, $\log p_{G}(\mathtt{tr}_{1})=\log p_{G}(\mathtt{tr}_{2})$ and $\log p_{\alpha}(\mathtt{tr}_{1})=\log p_{\alpha}(\mathtt{tr}_{2})=0$. If there is such a sample node $M_{\alpha}$, by assumption it is the only one whose address expression evaluations to $\alpha$ and by Theorem 1, $\log p_{G}(\mathtt{tr}_{1})-\log p_{G}(\mathtt{tr}_{2})=\log\bar{p}_{G_{k}}(\mathtt{tr}_{1},\sigma)-\log\bar{p}_{G_{k}}(\mathtt{tr}_{2},\sigma)$. After rearranging, we have

[TABLE]

∎

D.4. Proof of Corollary 3 - Correctness of LMH Optimisation

To prove correctness of the LMH optimisation, we require the notion of keys and minimal traces.

Definition 0.

The set of keys of a trace $\mathtt{tr}$ is the set of addresses that map to an non-$\mathtt{null}$ value

[TABLE]

A trace $\mathtt{tr}$ is minimal with respect to $G$ if

[TABLE]

Note that the density is undefined if a trace is read at an address that maps to $\mathtt{null}$, see (2). Since a trace is minimal if setting a single additional address to $\mathtt{null}$ results in undefined density, the minimality condition captures ”key errors”. For a minimal trace $\textnormal{keys}(\mathtt{tr})$ is the minimal set of addresses that need to be set to non-$\mathtt{null}$ values such that the density is well-defined, i.e. the trace contains no redundant information. We can immediately establish following lemma. Here, and in the following, we often omit the dependency on $G$ and write $p(\mathtt{tr})=p_{G}(\mathtt{tr})$ to declutter notation.

Lemma 0.

Let $\mathtt{tr}_{1}$ be a minimal trace and $\mathtt{tr}_{2}$ a trace such that $\forall\alpha\in\textnormal{keys}(\mathtt{tr}_{1})\colon\mathtt{tr}_{1}(\alpha)=\mathtt{tr}_{2}(\alpha)$. Then, the execution sequences of $\mathtt{tr}_{1}$ and $\mathtt{tr}_{2}$ as in Proposition 4 are the same and $p(\mathtt{tr}_{1})=p(\mathtt{tr}_{2}).$

Proof.

By induction on the execution sequence. Assume the execution sequences are equal up to sample node $N$ with address expression $E_{0}$ and state $\sigma$. Since $p(\mathtt{tr}_{1})\neq\textnormal{undefined}$, we must have $\mathtt{tr}_{1}(\sigma(E_{0}))\neq\mathtt{null}$ in (2). Thus, we have $\sigma(E_{0})\in\textnormal{keys}(\mathtt{tr}_{1})$ and $\mathtt{tr}_{2}(\sigma(E_{0}))=\mathtt{tr}_{1}(\sigma(E_{0}))\neq\mathtt{null}$ and the execution sequence continue to be equal. All other transition rules follow immediately. ∎

Next, we formally define the semantics of the LMH kernel (Wingate et al., 2011) in our setting.

Definition 0.

Let $Q(\beta)$ be a proposal distribution for each address $\beta$ and $\texttt{vs}(\beta)\sim Q(\beta)$ a sample from $Q(\beta)$. For a minimal trace $\mathtt{tr}$, the semantics of the LMH kernel (Wingate et al., 2011) that resamples at address $\alpha$ is given below, where we introduce new reserved variables $\boldsymbol{tr}$ and $\boldsymbol{q}$.

[TABLE]

Let there exists a state $\sigma\in\Sigma$ such that $(\sigma_{0},S)\Downarrow^{\mathtt{tr},\texttt{vs}}_{\texttt{LMH}(\alpha)}\sigma$ for program $S$ and initial program state $\sigma_{0}=(x_{1}\mapsto\mathtt{null},\dots,x_{n}\mapsto\mathtt{null},\textbf{p}\mapsto 1,\boldsymbol{tr}\mapsto\texttt{nulltr},\boldsymbol{q}\mapsto 1)$, where nulltr is a trace that maps all addresses to $\mathtt{null}$. Then, the new trace $\mathtt{tr}^{\prime}\coloneqq\sigma(\boldsymbol{tr})$ is called valid LMH proposal and $\mathtt{tr}^{\prime}$ is a minimal trace that agrees with $\mathtt{tr}$ on $\textnormal{keys}(\mathtt{tr})\cap\textnormal{keys}(\mathtt{tr}^{\prime})\setminus\{\alpha\}$. Furthermore, there exists a state $\sigma^{\prime}$ such that $(\sigma_{0},S)\Downarrow^{\mathtt{tr}^{\prime},\mathtt{tr}}_{\texttt{LMH}(\alpha)}\sigma^{\prime}$ and $\mathtt{tr}=\sigma^{\prime}(\boldsymbol{tr})$.

With $p(\mathtt{tr}^{\prime})=\sigma(\textbf{p})$, $Q_{\alpha}(\mathtt{tr}^{\prime}|\mathtt{tr})=\sigma(\boldsymbol{q})$, $p(\mathtt{tr})=\sigma^{\prime}(\textbf{p})$, and $Q_{\alpha}(\mathtt{tr}|\mathtt{tr}^{\prime})=\sigma^{\prime}(\boldsymbol{q})$ the acceptance probability is given by

[TABLE]

Lastly, let the equivalent semantics on the CFG be denoted by $\underset{\texttt{LMH}(\alpha)}{\overset{\mathtt{tr},\texttt{vs}}{\longrightarrow}}$.

Definition 0.

We can translate the LMH semantics to the semantics of context ForwardCtx, denoted by $\Downarrow^{\mathtt{tr},\texttt{vs}}_{\texttt{Fwd}(\alpha)}$, and of context BackwardCtx, denoted by $\Downarrow^{\mathtt{tr}}_{\texttt{Bck}(\alpha)}$.

Using log-densities, at the visit node, we have:

[TABLE]

At read nodes, we just read the trace $\mathtt{tr}$:

[TABLE]

The semantics of score nodes look similar to the semantics of sample statements in Definition 7 with the difference that $V_{0}\neq\alpha$:

[TABLE]

The semantics of the BackwardCtx are simply

[TABLE]

Now, we are ready to prove Corollary 3 which we make more precise here.

Statement.* For address $\alpha$ and minimal trace $\mathtt{tr}$, let vs be values from the proposal distributions as in Definition 7 and let $\mathtt{tr}^{\prime}$ be a valid LMH proposal. Assume that the address expression of at most one sample node in the execution sequence evaluates to $\alpha$ for all traces. Then, the execution sequences of $\mathtt{tr}$ and $\mathtt{tr}^{\prime}$ are equal up to $(\sigma,M_{\alpha})$, where $M_{\alpha}$ is the first and only sample node whose address expression evaluates to $\alpha$. For the sliced CFG $G_{k}$, where $\iota_{\textnormal{cfg}}(M_{\alpha})=N_{k}$, we have*

[TABLE]

*The quantities $\log\bar{p}_{G_{k}}(\mathtt{tr}^{\prime},\sigma)$, and $\log Q_{\alpha}(\mathtt{tr}^{\prime}|\mathtt{tr})$ are computed by ForwardCtx and the quantities $\log\bar{p}_{G_{k}}(\mathtt{tr},\sigma)$ and $\log Q_{\alpha}(\mathtt{tr}|\mathtt{tr}^{\prime})$ are computed by BackwardCtx. The new trace $\mathtt{tr}^{\prime}$ can be constructed from $\mathtt{tr}$ and the output of both contexts. *

Proof.

Since $\mathtt{tr}^{\prime}$ is a valid LMH proposal for $\mathtt{tr}$, we have that both traces are minimal and $\forall\beta\in\textnormal{keys}(\mathtt{tr})\cap\textnormal{keys}(\mathtt{tr}^{\prime})\setminus\{\alpha\}\colon\mathtt{tr}(\beta)=\mathtt{tr}^{\prime}(\beta)$.

Let $\mathtt{tr}_{1}=\mathtt{tr}$ and $\mathtt{tr}_{2}=\mathtt{tr}^{\prime}$ and define two new traces

[TABLE]

It holds that $\forall\beta\in\textnormal{keys}(\mathtt{tr}_{i})\colon\overline{\mathtt{tr}_{i}}(\beta)=\mathtt{tr}_{i}(\beta)$ (Note that $\alpha\in\textnormal{keys}(\mathtt{tr}_{i})$). By Lemma 6, $\overline{\mathtt{tr}_{1}}$ has the same execution sequence as $\mathtt{tr}$ and $\overline{\mathtt{tr}_{2}}$ the same as $\mathtt{tr}^{\prime}$.

Since $\overline{\mathtt{tr}_{1}}$ differs from $\overline{\mathtt{tr}_{2}}$ only at $\alpha$ the by Theorem 1 we have

[TABLE]

As the LMH semantics do not change the computation of $\sigma(\textbf{p})$, we immediately conclude that the ForwardCtx computes $\bar{p}_{G_{k}}(\mathtt{tr}^{\prime},\sigma)$ and the BackwardCtx computes $\bar{p}_{G_{k}}(\mathtt{tr},\sigma)$.

As in the proof of Theorem 1, but with LMH semantics, let $M^{\prime}_{1}$ be the last sample node in the execution sequence of $\mathtt{tr}$ such that $\iota_{\textnormal{cfg}}(M^{\prime}_{1})\in G_{k}$ and likewise $M^{\prime}_{2}$ for $\mathtt{tr}^{\prime}$.

[TABLE]

[TABLE]

For sample nodes that come before $M_{\alpha}$ or after $M^{\prime}_{1}$ or $M^{\prime}_{2}$ the address expressions evaluate to some $\beta\in\textnormal{keys}(\mathtt{tr})\cap\textnormal{keys}(\mathtt{tr}^{\prime})\setminus\{\alpha\}$. Thus, they do not update $\boldsymbol{q}$ and $\log Q_{\alpha}(\mathtt{tr}^{\prime}|\mathtt{tr})=\log\sigma_{1}(\boldsymbol{q})$ as well as $\log Q_{\alpha}(\mathtt{tr}|\mathtt{tr}^{\prime})=\log\sigma_{2}(\boldsymbol{q})$ are fully computed in the sliced CFG by ForwardCtx and BackwardCtx, respectively (compare Definition 7).

For the same reason, the difference between $\mathtt{tr}$ and $\mathtt{tr}^{\prime}$ is also fully accounted for in the sliced CFG. The new trace $\mathtt{tr}^{\prime}$ can be constructed by removing the address-value pairs encountered in the BackwardCtx from $\mathtt{tr}$ and adding the address-value pairs encountered in the ForwardCtx.

∎

D.5. Proof of Corollary 4 - Correctness of BBVI optimisation

Before we can prove correctness of the BBVI optimisation, we need to formally define the variational distribution over traces $Q_{\phi}$.

Definition 0.

For every address $\alpha$, let there be a probability measure on some measurable space $\mathcal{M_{\alpha}}=(Z_{\alpha},\Sigma_{\alpha})$ parameterised by $\phi_{\alpha}$ that admits a density $q_{\alpha}(z_{\alpha}|\phi_{\alpha})$ differentiable with respect to $\phi_{\alpha}$ (the density is with respect to the Lebesgue measure on $\mathbb{R}^{d_{\alpha}}$ or counting measure on $\mathbb{Z}^{d_{\alpha}}$).

For every finite address set $A$, define the joint density $q_{A}((z_{\alpha})_{\alpha\in A}|(\phi_{\alpha})_{\alpha\in A})=\prod_{\alpha\in A}q_{\alpha}(z_{\alpha}|\phi_{\alpha})$ on the product space $\mathcal{M}_{A}=(Z_{A},\Sigma_{A})=\bigotimes_{\alpha\in A}\mathcal{M_{\alpha}}.$

On the set of traces $\mathcal{T}$, define the $\sigma$-algebra $\Sigma_{\mathcal{T}}$, such that $T\in\Sigma_{\mathcal{T}}$ if and only if for all finite address sets $A$, we have $\{(\mathtt{tr}(\alpha))_{\alpha\in A}\colon\mathtt{tr}\in T\}\in\Sigma_{A}$. Lastly, assume that for program density $p$ the set $\mathcal{T}_{\textnormal{minimal}}=\{\mathtt{tr}\in\mathcal{T}\colon\mathtt{tr}\text{ is minimal wrt. }p\}$ is measurable in $\Sigma_{\mathcal{T}}$.

Then, the variational distribution is defined by the density

[TABLE]

where $\mathcal{T}_{A}\coloneqq\{\mathtt{tr}\in\mathcal{T}\colon\forall\alpha\in\mathbf{Strings}\colon\mathtt{tr}(\alpha)\neq\mathtt{null}\Leftrightarrow\alpha\in A\}$. For all measurable functions $f$ we can express the expectation with respect to $Q_{\phi}$ as

[TABLE]

Above $\iota_{\mathcal{T}}(z_{A})$ for $z_{A}=(z_{\alpha})_{\alpha\in A}\in Z_{A}$ is the trace that maps address $\alpha$ to $z_{\alpha}$ if $\alpha\in A$ and all other addresses to $\mathtt{null}$. Note that we may exchange integration and sum, since for a minimal trace only one term in the sum is non-zero.

Before repeating the correctness statement in terms of Definition 9, we formally defined the semantics of BBVICtx.

Definition 0.

The semantics of the BBVICtx are very similar to the standard semantics with the only difference that we store $\log q_{\alpha}(\mathtt{tr}(\alpha)|\phi_{\alpha})$ at visit statements and the use of log-densities.

[TABLE]

The semantics of read statements remains unchanged:

[TABLE]

Lastly, the semantics of score statements are only rewritten in terms of log-densities:

[TABLE]

Statement.* Let $Q_{\phi}$ be a variational distribution over traces as in Definition 9. Assume that the address expression of at most one sample node in the execution sequence evaluates to $\alpha$ for all traces. Assume that the program density $p$ and the set $\mathcal{T}_{\textnormal{minimal}}$ are measurable with respect to $(\mathcal{T},\Sigma_{\mathcal{T}})$. With $p_{\alpha}$ as in Corollary 2, it holds that*

[TABLE]

*The quantity $\log p_{\alpha}(\mathtt{tr})-\log q_{\alpha}(\mathtt{tr}(\alpha)|\phi_{\alpha})$ is computed by BBVICtx. *

Proof.

First, we repeat the known identity for the BBVI estimator

[TABLE]

Then, for a minimal trace $\mathtt{tr}$, we have

[TABLE]

Further, let $\Delta_{p}(\mathtt{tr})=\log p(\mathtt{tr})-\log p_{\alpha}(\mathtt{tr})$ and $\Delta_{q}(\mathtt{tr})=\log Q_{\phi}(\mathtt{tr})-\log q_{\alpha}(\mathtt{tr}(\alpha)|\phi_{\alpha})$. For an address set $A$ with $\alpha\in A$, by Corollary 2 and by definition of $Q_{\phi}$, both $\Delta_{p}\circ\iota_{\mathcal{T}}$ and $\Delta_{q}\circ\iota_{\mathcal{T}}$ are constant with respect to $z_{\alpha}$ in $Z_{A}$. Thus,

[TABLE]

because \int\nabla_{\phi_{\alpha}}\big{[}\log q_{\alpha}(z_{\alpha}|\phi_{\alpha})\big{]}q_{\alpha}(z_{\alpha}|\phi_{\alpha})\,dz_{\alpha}=0. Similarly,

[TABLE]

Lastly, we rewrite

[TABLE]

Since we have shown that the second term in the sum integrates to [math] for every address set $A$, with $\alpha\in A$, by linearity of the integral, we have

[TABLE]

Finally, from the definition of the BBVICtx it follows immediately that this context computes the quantities $\log p_{\alpha}(\mathtt{tr})$ and $\log q_{\alpha}(\mathtt{tr}(\alpha)|\phi_{\alpha})$. ∎

D.6. Proof of Proposition 5 - Correctness of SMC optimisation

First, we formalise the sliced CFG construction for SMC.

Definition 0.

Let $G=(N_{\textnormal{start}},N_{\textnormal{end}},\mathcal{N},\mathcal{E})$ be a CFG with $K$ sample nodes, which we enumerate with ${N_{j}=\textnormal{Assign}(x_{j}=\mbox{\leavevmode\lstinline{{\lst@@@set@language\lst@@@set@numbers\lst@@@set@frame\lst@@@set@rulecolor{\@listingGroup{ltx_lst_identifier}{sample}}}}}}(E_{0}^{j},f^{j}(E_{1}^{j},\dots,E_{n_{j}}^{j})))$ for $j=1,\dots,K$.

We define the set of nodes that lie on path from $N_{k}$ to some of other sample node $N_{j}$:

[TABLE]

In this construction, we do not require read statements. We restrict the set of edges to the subset of nodes $\mathcal{N}_{k}$:

[TABLE]

We connect the new start node $N^{k}_{\textnormal{start}}$ to $N_{k}$ and the new end node $N^{k}_{\textnormal{end}}$ to all $N_{j}$, $j\in\{1,\dots,K\}$:

[TABLE]

Finally, define the CFG sliced with respect to $N_{k}$ as

[TABLE]

The ”continuation property” of this sliced CFG follows directly from its construction.

Statement.* Let $G$ be a CFG such that for initial program state $\sigma_{0}$ and trace $\mathtt{tr}$, we have execution sequence $(\sigma_{0},\text{START})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma_{i},M_{i})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma_{l},\text{END})$ in its unrolled version $H$. Let $(\sigma^{\prime}_{1},M^{\prime}_{1})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{2},M^{\prime}_{2})$ be a sub-sequence such that $\iota_{\textnormal{cfg}}(M^{\prime}_{1})=N_{k}$ and $\iota_{\textnormal{cfg}}(M^{\prime}_{2})=N_{j}$ are the only sample nodes. Then, the execution sequence in the unrolled version of the sliced CFG $G_{k}^{\texttt{SMC}}$ matches the sub-sequence: $(\sigma^{\prime}_{1},\text{START})\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{1},M^{\prime}_{1})\overset{\mathtt{tr}}{\longrightarrow}\cdots\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{2},M^{\prime}_{2})\overset{\mathtt{tr}}{\longrightarrow}(\sigma^{\prime}_{2},\text{END}).$ *

Proof.

Follows directly from Definition 11. ∎