Towards Verified Stochastic Variational Inference for Probabilistic Programs

TL;DR

This paper analyzes the score estimator variational inference algorithm used in probabilistic programming, formalizes its assumptions, and develops static analysis tools to verify these assumptions automatically, improving reliability in deep probabilistic languages like Pyro.

Contribution

It provides a formal semantics for the score estimator, identifies implicit assumptions, and introduces automated static analysis rules to verify these assumptions in probabilistic programs.

Findings

Identified a bug in one Pyro model-guide pair

Revealed non-standard inference engine use in another case

Verified assumptions in remaining cases

Abstract

Probabilistic programming is the idea of writing models from statistics and machine learning using program notations and reasoning about these models using generic inference engines. Recently its combination with deep learning has been explored intensely, leading to the development of deep probabilistic programming languages such as Pyro. At the core of this development lie inference engines based on stochastic variational inference algorithms. When asked to find information about the posterior distribution of a model written in such a language, these algorithms convert this posterior-inference query into an optimisation problem and solve it approximately by gradient ascent. In this paper, we analyse one of the most fundamental and versatile variational inference algorithms, called score estimator, using tools from denotational semantics and program analysis. We formally express what this algorithm does on models denoted by programs, and expose implicit assumptions made by the algorithm. The violation of these assumptions may lead to an undefined optimisation objective or the loss of convergence guarantee of the optimisation process. We then describe rules for proving these assumptions, which can be automated by static program analyses. Some of our rules use nontrivial facts from continuous mathematics, and let us replace requirements about integrals in the assumptions, by conditions involving differentiation or boundedness, which are much easier to prove automatically. Following our general methodology, we have developed a static program analysis for Pyro that aims at discharging the assumption about what we call model-guide support match. Applied to the eight representative model-guide pairs from the Pyro webpage, our analysis finds a bug in one of these cases, reveals a non-standard use of an inference engine in another, and shows the assumptions are met in the remaining cases.

Figures

Tables

Table 1. Table 1. Key features of the model-guide pairs from Pyro examples. LoC denotes the lines of code of model and guide. The columns “Total #” show the number of objects/commands of each type used in model and guide, and the columns “Total dimension” show the total dimension of tensors in model and guide, either sampled from sample or used inside score , as well as the dimension of θ 𝜃 \theta in guide.

			Total #				Total dimension
Name	Corresponding probabilistic model	LoC	for	plate	sample	score	sample	score	$\theta$
br	Bayesian regression	27	0	1	10	1	10	170	9
csis	Inference compilation	31	0	0	2	2	2	2	480
lda	Amortised latent Dirichlet allocation	76	0	5	8	1	21008	64000	121400
vae	Variational autoencoder (VAE)	91	0	2	2	1	25600	200704	353600
sgdef	Deep exponential family	94	0	8	12	1	231280	1310720	231280
dmm	Deep Markov model	246	3	2	2	1	640000	281600	594000
ssvae	Semi-supervised VAE	349	0	2	4	1	24000	156800	844000
air	Attend-infer-repeat	410	2	2	6	1	20736	160000	6040859

Table 2. (a) Results for Pyro test suite. 39 model-guide pairs are grouped into 7 categories, based on which type of plate objects are used. #Same (or #Diff) denotes the number of model-guide pairs for which the output of our analyser, valid or invalid, is the same as (or different from) the documented output. #Crash denotes the number of pairs for which our analyser crashes.

Category	#Same	#Diff	#Crash	Time
No plates	9	0	0	0.001
Single for-plate	4	0	3	0.004
Nested for-plates	2	0	2	0.026
Single with-plate	5	0	0	0.001
Nested with-plates	7	2	0	0.002
Non-nested with-plates	2	0	0	0.002
Nested for-plate & with-plate	0	0	3	N/A
Total	29	2	8	0.003

Table 3. (a) Results for Pyro test suite. 39 model-guide pairs are grouped into 7 categories, based on which type of plate objects are used. #Same (or #Diff) denotes the number of model-guide pairs for which the output of our analyser, valid or invalid, is the same as (or different from) the documented output. #Crash denotes the number of pairs for which our analyser crashes.

Category	#Same	#Diff	#Crash	Time
No plates	9	0	0	0.001
Single for-plate	4	0	3	0.004
Nested for-plates	2	0	2	0.026
Single with-plate	5	0	0	0.001
Nested with-plates	7	2	0	0.002
Non-nested with-plates	2	0	0	0.002
Nested for-plate & with-plate	0	0	3	N/A
Total	29	2	8	0.003

Table 4. (b) Results for Pyro examples. The column “Valid?” shows the output of our analysis, valid or invalid.

Name	Valid?	Time
br	x	0.006
csis	o	0.007
lda	x	0.014
vae	o	0.005
sgdef	o	0.070
dmm	o	0.536
ssvae	o	0.013
air	o	4.093

Equations

$$p(v|obs=0)=\frac{p(v,obs=0)}{\int \mathrm{d}vp(v,obs=0)}=\frac{p(v)\cdot p(obs=0|v)}{p(obs=0)}.$$

$${\mathrm{argmin}}_{\theta }\mathrm{KL}(q_{\theta }(v)||p(v|obs=0)),\text{where }\mathrm{KL}(q_{\theta }(v)||p(v|obs=0))\triangleq {\mathbb{E}}_{q_{\theta }(v)}\left[\log \frac{q_{\theta }(v)}{p(v|obs=0)}\right].$$

$$\nabla_{\theta}\mathrm{KL}(q_{\theta}(v)||p(v|\mathit{obs}{=}0))\Big{|}_{\theta=\theta_{n}}\approx\frac{1}{N}\sum_{i=1}^{N}(\nabla_{\theta}\log q_{\theta_{n}}(v_{i}))\cdot\log\frac{q_{\theta_{n}}(v_{i})}{p(v_{i},\mathit{obs}{=}0)}$$

$${\theta }_{n+1}\leftarrow {\theta }_n-0.01\times \frac{1}{N}\sum _{i=1}^N({\nabla }_{\theta }\log q_{{\theta }_n}(v_i))\cdot \log \frac{q_{{\theta }_n}(v_i)}{p(v_i,obs=0)}.$$

$$\nabla_{\theta}\mathrm{KL}(q_{\theta}(v)||p(v|\mathit{obs}{=}0))\Big{|}_{\theta=\theta_{n}}=\mathbb{E}\left[{\frac{1}{N}\sum_{i=1}^{N}(\nabla_{\theta}\log q_{\theta_{n}}(v_{i}))\cdot\log\frac{q_{\theta_{n}}(v_{i})}{p(v_{i},\mathit{obs}{=}0)}}\right]$$

$$\mu\Big{(}\bigcup_{n=0}^{\infty}B_{n}\Big{)}=\sum_{n=0}^{\infty}\mu(B_{n}).$$

$$\bigotimes_{i\in I}\Sigma_{i}\triangleq\sigma\Big{(}\Big{\{}\prod_{i}A_{i}\,\Big{|}~{}\,\text{ there is a finite $I_{0}\subseteq I$ such that }\ (\forall j\in I\setminus I_{0}.\,A_{j}=X_{j})\wedge(\forall i\in I_{0}.\,A_{i}\in\Sigma_{i})\Big{\}}\Big{)}.$$

$$\underset{i\in I}{⨁}{\Sigma }_i\triangleq \{A\mid A\cap X_i\in {\Sigma }_i\text{ for all i∈I}\}.$$

$$\mathrm{Mea}(X)$$

$$\Big{(}\bigotimes_{i\in I}\mu_{i}\Big{)}\Big{(}\prod_{i\in I}A_{i}\Big{)}=\prod_{i\in I}\mu_{i}(A_{i}).$$

$$\Big{(}\bigoplus_{i\in I}\mu_{i}\Big{)}\Big{(}\bigcup_{i\in I}A_{i}\Big{)}=\sum_{i\in I}\mu_{i}(A_{i}).$$

$$\begin{array}[]{@{}r@{}l@{\qquad\quad}rl@{}}\text{\it Real constants}\quad c&\;\;\in\;\mathbb{R}&\text{\it Primitive real-valued functions}\quad f\;::=\;\ldots\\[2.15277pt] \text{\it String constants}\quad\alpha&\;\;\in\;\mathit{{Str}}&\text{\it Primitive string-valued functions}\quad g\;::=\;\ldots\\[2.15277pt] \text{\it Real expressions}\quad E&\hbox to0.0pt{${}\;::=\;c\,\mid\,x\,\mid\,f(E,\ldots,E)$\hss}\\[2.15277pt] \text{\it Boolean expressions}\quad B&\hbox to0.0pt{${}\;::=\;\mathsf{{true}}\,\mid\,E<E\,\mid\,B\wedge B\,\mid\,\neg B$\hss}\\[2.15277pt] \text{\it String expressions}\quad S&\hbox to0.0pt{${}\;::=\;\alpha\,\mid\,g(S,\ldots,S,E,\ldots,E)$\hss}\\[2.15277pt] \text{\it Commands}\quad C&\hbox to0.0pt{${}\;::=\;\mathsf{{skip}}\,\mid\,x:=E\,\mid\,C;C\,\mid\,\mathsf{{if}}\,B\,\{C\}\,\mathsf{{else}}\,\{C\}\,\mid\,\mathsf{{while}}\,B\,\{C\}$\hss}\\[2.15277pt] &\hbox to0.0pt{${}~{}\;\;\mid\;\;~{}x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E,E)\,\mid\,\mathsf{{score}}_{\mathsf{{norm}}}(E,E,E)$\hss}\end{array}$$

$$x:={\mathsf{sample}}_{\mathsf{norm}}(‘‘a",0.0,5.0);{\mathsf{score}}_{\mathsf{norm}}(3.0,x,1.0).$$

$$[[E]]$$

$$\text{Stores}s$$

$$\text{Random databases}r$$

$$\text{States}\sigma$$

$$[[x]]s$$

$$[[\alpha ]]s$$

$$\begin{array}[]{@{}r@{\;}c@{\;}l@{}}{\llbracket\mathsf{{skip}}\rrbracket}(\sigma)(A)&\triangleq&{\mathds{1}_{[{(\sigma,1)\in A}]}}\\[4.30554pt] {\llbracket x:=E\rrbracket}(\sigma)(A)&\triangleq&{\mathds{1}_{[{(\sigma[x\mapsto{\llbracket E\rrbracket}\sigma_{s}],1)\in A}]}}\\[4.30554pt] {\llbracket C_{0};C_{1}\rrbracket}(\sigma)(A)&\triangleq&\int{\llbracket C_{0}\rrbracket}(\sigma)(\mathrm{d}(\sigma^{\prime},w^{\prime}))\int{\llbracket C_{1}\rrbracket}(\sigma^{\prime})(\mathrm{d}(\sigma^{\prime\prime},w^{\prime\prime}))\,{\mathds{1}_{[{(\sigma^{\prime\prime},w^{\prime}\cdot w^{\prime\prime})\in A}]}}\\[4.30554pt] {\llbracket\mathsf{{if}}\,B\,\{C_{0}\}\,\mathsf{{else}}\,\{C_{1}\}\rrbracket}(\sigma)(A)&\triangleq&{\mathds{1}_{[{{\llbracket B\rrbracket}\sigma_{s}=\mathit{{true}}}]}}\cdot{\llbracket C_{0}\rrbracket}(\sigma)(A)+{\mathds{1}_{[{{\llbracket B\rrbracket}\sigma_{s}\neq\mathit{{true}}}]}}\cdot{\llbracket C_{1}\rrbracket}(\sigma)(A)\\[4.30554pt] {\llbracket\mathsf{{while}}\,B\,\{C\}\rrbracket}(\sigma)(A)&\triangleq&(\mathit{{fix}}\,F)(\sigma)(A)\\ \lx@intercol\qquad\big{(}\text{where }F(\kappa)(\sigma)(A)\triangleq{\mathds{1}_{[{{\llbracket B\rrbracket}\sigma_{s}\neq\mathit{{true}}}]}}\cdot{\mathds{1}_{[{(\sigma,1)\in A}]}}\hfil\\ \lx@intercol\phantom{\qquad\big{(}\text{where }F(\kappa)(\sigma)(A)\triangleq}{}+{\mathds{1}_{[{{\llbracket B\rrbracket}\sigma_{s}=\mathit{{true}}}]}}\cdot\int{\llbracket C\rrbracket}(\sigma)(\mathrm{d}(\sigma^{\prime},w^{\prime}))\int\kappa(\sigma^{\prime})(\mathrm{d}(\sigma^{\prime\prime},w^{\prime\prime}))\,{\mathds{1}_{[{(\sigma^{\prime\prime},w^{\prime}\cdot w^{\prime\prime})\in A}]}}\big{)}\hfil\\[4.30554pt] {\llbracket x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})\rrbracket}(\sigma)(A)&\triangleq&{}{\mathds{1}_{[{{\llbracket S\rrbracket}\sigma_{s}\not\in\mathrm{dom}(\sigma_{r})}]}}\cdot{\mathds{1}_{[{{\llbracket E_{2}\rrbracket}\sigma_{s}\in(0,\infty)}]}}\\[2.15277pt] &&\cdot\int\mathrm{d}v\,\big{(}\mathcal{N}(v;{\llbracket E_{1}\rrbracket}\sigma_{s},{\llbracket E_{2}\rrbracket}\sigma_{s})\cdot{\mathds{1}_{[{((\sigma_{s}[x\mapsto v],\sigma_{r}[{\llbracket S\rrbracket}\sigma_{s}\mapsto v]),1)\in A}]}}\big{)}\\[4.30554pt] {\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}(\sigma)(A)&\triangleq&{\mathds{1}_{[{{\llbracket E_{2}\rrbracket}\sigma_{s}\in(0,\infty)}]}}\cdot{\mathds{1}_{[{(\sigma,\mathcal{N}({\llbracket E_{0}\rrbracket}\sigma_{s};{\llbracket E_{1}\rrbracket}\sigma_{s},{\llbracket E_{2}\rrbracket}\sigma_{s}))\in A}]}}\end{array}$$

$$\mathrm{Mea}(C,A)$$

$$r\#r^{\prime }⟺\mathrm{dom}(r)\cap \mathrm{dom}(r^{\prime })=\varnothing .$$

$$\mathrm{dom}(r\uplus r^{\prime })\triangleq \mathrm{dom}(r)\cup \mathrm{dom}(r^{\prime });(r\uplus r^{\prime })(\alpha )\triangleq \text{if }\alpha \in \mathrm{dom}(r)\text{ then }r(\alpha )\text{ else }{r}^{\prime }(\alpha ).$$

$$\rho(R)\triangleq\sum_{\begin{subarray}{c}K\subseteq_{\mathit{{fin}}}\mathit{{Str}}\end{subarray}}\Big{(}\bigotimes_{\alpha\in K}\lambda\Big{)}(R\cap[K\to\mathbb{R}]),$$

$$g(s,r)=(s^{\prime },r^{\prime },w^{\prime },p^{\prime })⟹$$

$$\wedge (\forall r^{\prime \prime \prime }.r\#r^{\prime \prime \prime }⟹g(s,r\uplus r^{\prime \prime \prime })=(s^{\prime },r^{\prime }\uplus r^{\prime \prime \prime },w^{\prime },p^{\prime })).$$

$$\mathcal{D}\triangleq$$

$$g⊑g^{\prime }⟺\forall (s,r)\in Store\times RDB.(g(s,r)=\perp \vee g(s,r)=g^{\prime }(s,r)).$$

$$\displaystyle g^{\ddagger}(\bot)\triangleq\bot,\quad\;g^{\ddagger}(s,r,w,p)\triangleq\left\{\begin{array}[]{ll}\bot&\text{ if }g(s,r)=\bot,\\[2.15277pt] (s^{\prime},r^{\prime},w\times w^{\prime},p\times p^{\prime})&\text{ if }g(s,r)=(s^{\prime},r^{\prime},w^{\prime},p^{\prime}).\end{array}\right.$$

$$\begin{array}[]{@{}r@{\;}c@{\;}l@{}}{\llbracket\mathsf{{skip}}\rrbracket}_{d}(s,r)&\triangleq&(s,r,1,1)\\[4.30554pt] {\llbracket x{:=}E\rrbracket}_{d}(s,r)&\triangleq&(s[x\mapsto{\llbracket E\rrbracket}s],r,1,1)\\[4.30554pt] {\llbracket C_{0};C_{1}\rrbracket}_{d}(s,r)&\triangleq&({\llbracket C_{1}\rrbracket}_{d}^{\ddagger}\circ{\llbracket C_{0}\rrbracket}_{d})(s,r)\\[4.30554pt] {\llbracket\mathsf{{if}}\,B\,\{C_{0}\}\,\mathsf{{else}}\,\{C_{1}\}\rrbracket}_{d}(s,r)&\triangleq&\text{if}\ ({\llbracket B\rrbracket}s=\mathit{{true}})\ \text{then}\ {\llbracket C_{0}\rrbracket}_{d}(s,r)\ \text{else}\ {\llbracket C_{1}\rrbracket}_{d}(s,r)\\[4.30554pt] {\llbracket\mathsf{{while}}\,B\,\{C\}\rrbracket}_{d}(s,r)&\triangleq&(\mathit{{fix}}\,G)(s,r)\\ \lx@intercol\qquad\qquad\qquad\qquad\qquad(\text{where }G(g)(s,r)=\text{if}\;({\llbracket B\rrbracket}s\neq\mathit{{true}})\;\text{then}\;(s,r,1,1)\;\text{else}\;(g^{\ddagger}\circ{\llbracket C\rrbracket}_{d})(s,r))\hfil\\[4.30554pt] {\llbracket x{:=}\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})\rrbracket}_{d}(s,r)&\triangleq&\text{if}\;({\llbracket S\rrbracket}s\not\in\mathrm{dom}(r)\vee{\llbracket E_{2}\rrbracket}s\not\in(0,\infty))\;\text{then}\;\bot\\ &&\text{else}\;(s[x\mapsto r({\llbracket S\rrbracket}s)],\,r\setminus{\llbracket S\rrbracket}s,\,1,\,\mathcal{N}(r({\llbracket S\rrbracket}s);{\llbracket E_{1}\rrbracket}s,{\llbracket E_{2}\rrbracket}s))\\[4.30554pt] {\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}_{d}(s,r)&\triangleq&\text{if}\;({\llbracket E_{2}\rrbracket}s\not\in(0,\infty))\;\text{then}\;\bot\;\text{else}\;(s,r,\mathcal{N}({\llbracket E_{0}\rrbracket}s;{\llbracket E_{1}\rrbracket}s,{\llbracket E_{2}\rrbracket}s),1)\end{array}$$

$$dens(g,s)$$

$$get(g,s)$$

$$\displaystyle\int\rho(\mathrm{d}r)\,\Big{(}\mathit{{dens}}(g_{2}^{\ddagger}\circ g_{1},s)(r)\cdot h\Big{(}\mathit{{get}}(g_{2}^{\ddagger}\circ g_{1},s)(r),r\Big{)}\Big{)}$$

Full text

probabilistic programming, static analysis, semantics and correctness

††journal: PACMPL††journalvolume: 0††journalnumber: POPL††article: 0††journalyear: 2020††publicationmonth: 1††doi: 10.1145/nnnnnnn.nnnnnnn††copyright: none††ccs: Mathematics of computing Bayesian computation††ccs: Mathematics of computing Variational methods††ccs: Theory of computation Probabilistic computation††ccs: Theory of computation Denotational semantics††ccs: Software and its engineering Correctness††ccs: Software and its engineering Automated static analysis††ccs: Computing methodologies Machine learning

1. Introduction

Probabilistic programming refers to the idea of writing models from statistics and machine learning using program notations and reasoning about these models using generic inference engines. It has been the subject of active research in machine learning and programming languages, because of its potential for enabling scientists and engineers to design and explore sophisticated models easily; when using these languages, they no longer have to worry about developing custom inference engines for their models, a highly-nontrivial task requiring expertise in statistics and machine learning. Several practical probabilistic programming languages now exist, and are used for a wide range of applications (Carpenter et al., 2017; Minka et al., 2014; Gordon et al., 2014; Goodman et al., 2008; Mansinghka et al., 2014; Wood et al., 2014; Tolpin et al., 2016; Narayanan et al., 2016; Gehr et al., 2016).

In this paper, we consider inference engines that lie at the core of so called deep probabilistic programming languages, such as Pyro (Bingham et al., 2019), Edward (Tran et al., 2016, 2018) and ProbTorch (Siddharth et al., 2017). These languages let their users freely mix deep neural networks with constructs from probabilistic programming, in particular, those for writing Bayesian probabilistic models. In so doing, they facilitate the development of probabilistic deep-network models that may address the problem of measuring the uncertainty in current non-Bayesian deep-network models; a non-Bayesian model may predict that the price of energy goes up and that of a house goes down, but it cannot express, for instance, that the model is very confident with the first prediction but not the second.

The primary inference engines for these deep probabilistic programming languages implement stochastic (or black-box) variational inference111The term stochastic variational inference (VI) often refers to VI with data subsampling (Hoffman et al., 2013), and our usage of the term is often called black-box VI (Ranganath et al., 2014) to stress the treatment of a model as a black-box sampler. algorithms. Converting inference problems into optimisation problems is the high-level idea of these algorithms.222The inference problems in their original forms involve solving summation/integration/counting problems, which are typically more difficult than optimisation problems. The variational-inference algorithms convert the former problems to the latter ones, by looking for approximate, not exact, answers to the former. When asked to find information about the posterior distribution of a model written in such a language, these algorithms convert this question to an optimisation problem and solve the problem approximately by performing a gradient descent or ascent on the optimisation objective. The algorithms work smoothly with gradient-based parameter-learning algorithms for deep neural networks, which is why they form the backbone for deep probabilistic programming languages.

In this paper, we analyse one of the most fundamental and versatile variational inference algorithms, called score estimator or REINFORCE333REINFORCE (Williams, 1992) is an algorithm originally developed for reinforcement learning (RL), but it is commonly used as a synonym of the score-estimator algorithm. This is because REINFORCE and score estimator use a nearly identical method for estimating the gradient of an optimisation objective. (Williams, 1992; Paisley et al., 2012; Wingate and Weber, 2013; Ranganath et al., 2014), using tools from denotational semantics and program analysis (Cousot and Cousot, 1977, 1979, 1992). We formally express what this algorithm does on models denoted by probabilistic programs, and expose implicit assumptions made by the algorithm on the models. The violation of these assumptions can lead to undefined optimisation objective or the loss of convergence guarantee of the optimisation process. We then describe rules for proving these assumptions, which can be automated by static program analyses. Some of our rules use nontrivial facts from continuous mathematics, and let us replace requirements about integrals in the assumptions, such as integrability of functions defined in terms of programs’ denotations, by the conditions involving differentiation or boundedness, which are much easier to prove automatically (and manually) than the original requirements.

Following our general methodology, we have developed a static program analysis for the Pyro programming language that can discharge one assumption of the inference algorithm about so called model-guide pairs. In Pyro and other deep probabilistic programming languages, a program denoting a model typically comes with a companion program, called guide, decoder, or inference network. This companion, which we call guide, helps the inference algorithm to find a good approximation to what the model ultimately denotes under a given dataset (i.e., the posterior distribution of the model under the dataset); the algorithm uses the guide to fix the search space of approximations, and solves an optimisation problem defined on that space. A model and a guide should satisfy an important correspondence property, which says that they should use the same sets of random variables, and for any such random variable, if the probability of the variable having a particular value is zero in the model, it should also be zero in the guide. If the property is violated, the inference algorithm may attempt to solve an optimisation problem with undefined optimisation objective and return parameter values that do not make any sense. Our static analysis checks this correspondence property for Pyro programs. When applied to eight representative model-guide pairs from the Pyro webpage, which include sophisticated neural network models such as Attend-Infer-Repeat (AIR), the analysis found a bug in one of these cases, revealed a non-standard use of the inference algorithm in another, and proved that the property holds in the remaining six cases.

Another motivation for this paper is to demonstrate an opportunity for programming languages and verification research to have an impact on the advances of machine learning and AI technologies. One popular question is: what properties should we verify on machine-learning programs? Multiple answers have been proposed, which led to excellent research results, such as those on robustness of neural networks (Mirman et al., 2018). But most of the existing research focuses on the final outcome of machine learning algorithms, not the process of applying these algorithms. One of our main objectives is to show that the process often relies on multiple assumptions on models and finding automatic ways for discharging these assumptions can be another way of making PL and verification techniques contribute. While our suggested solutions are not complete, they are intended to show the richness of this type of problems in terms of theory and practice.

We summarise the contributions of the paper:

• •

We formally express the behaviour of the most fundamental variational inference algorithm on probabilistic programs using denotational semantics, and identify requirements on program denotations that are needed for this algorithm to work correctly.

• •

We describe conditions that imply the identified requirements but are easier to prove. The sufficiency of the conditions relies on nontrivial results from continuous mathematics. We sketch a recipe for building program analyses for checking these conditions automatically.

• •

We present a static analysis for the Pyro language that checks the correspondence requirement of model-guide pairs. The analysis is based on our recipe, but extends it significantly to address challenges for dealing with features of the real-world language. Our analysis has successfully verified 6 representative Pyro model-guide examples, and found a bug in one example.

The proofs of lemmas and theorems of the paper are included in Appendix.

2. Variational Inference and Verification Challenges by Examples

We start by explaining informally the idea of stochastic variational inference (in short SVI), one fundamental SVI algorithm, and the verification challenges that arise when we use this algorithm.

2.1. Stochastic variational inference

In a probabilistic programming language, we specify a model by a program. The program model() in Figure 1(a) is an example. It describes a joint probability density $p(v,\mathit{obs})$ on two real-valued random variables $v$ and $\mathit{obs}$. The value of the former is not observed, while the latter is observed to have the value [math]. Finding out the value of $v$ is the objective of writing this model. The joint density $p(v,\mathit{obs})$ is expressed in terms of prior $p(v)$ and likelihood $p(\mathit{obs}|v)$ in the program. The prior $p(v)$ of $v$ is the normal distribution with mean [math] and standard deviation $5$, and it expresses the belief about the possible value of $v$ before any observation. The likelihood $p(\mathit{obs}|v)$ is a normal distribution whose mean and standard deviation are either $(1,1)$ or $(-2,1)$ depending on the sign of the value of $v$. The purpose of most inference algorithms is to compute exactly or approximately the posterior density given a prior and a likelihood. In our example, the posterior $p(v|\mathit{obs}{=}0)$ is:

[TABLE]

Intuitively, the posterior expresses an updated belief on $v$ upon observing $\mathit{obs}=0$. The dashed blue and solid orange lines in Figure 1(b) show the prior and posterior densities, respectively. Note that the density of a positive $v$ in the prior went up in the posterior. This is because when $v>0$, the mean of $p(\mathit{obs}|v)$ is $1$, a value closer to the observed value [math] than the alternative $-2$ for the mean.

SVI algorithms approach the posterior inference problem from the optimisation angle. They consider a collection of approximating distributions to a target posterior, formulate the problem of finding a good approximation in the collection as an optimisation problem, and solve the optimisation problem. The solution becomes the result of those algorithms. In Pyro, the users specify such a collection by a single parameterised program called guide; the collection can be generated by instantiating the parameters with different values. The program guide() in Figure 1(a) is an example. It has a real-valued parameter $\theta$ (written as theta in the program), and states that the probability density $q_{\theta}(v)$ of $v$ is the normal distribution with unknown mean $\theta$ and standard deviation $1$. The lines 13–17 in the figure show how to apply a standard SVI engine of Pyro (called Trace_ELBO) to find a good $\theta$. They instruct the engine to solve the following optimisation problem:

[TABLE]

The optimisation objective $\mathrm{KL}(q_{\theta}(v)||p(v|\mathit{obs}{=}0))$ is the KL divergence from $q_{\theta}(v)$ to $p(v|\mathit{obs}{=}0)$, and measures the similarity between the two densities, having a small value when the densities are similar. The KL divergence is drawn in Figure 1(c) as a function of $\theta$, and the dotted green line in Figure 1(b) draws the density $q_{\theta}$ at the optimum $\theta$. Note that the mean of this distribution is biased toward the positive side, which reflects the fact that the property $v>0$ has a higher probability than its negation $v\leq 0$ in the posterior distribution.

One of the most fundamental and versatile algorithms for SVI is score estimator (also called REINFORCE). It repeatedly improves $\theta$ in two steps. First, it estimates the gradient of the optimisation objective with samples from the current $q_{\theta_{n}}$:

[TABLE]

where $v_{1},\ldots,v_{N}$ are independent samples from the distribution $q_{\theta_{n}}$. Then, the algorithm updates $\theta$ with the estimated gradient (the specific learning rate $0.01$ is chosen to improve readability):

[TABLE]

When the learning rate $0.01$ is adjusted according to a known scheme, the algorithm is guaranteed to converge to a local optimum (in many cases) because its gradient estimate satisfies the following unbiasedness property (in those cases):

[TABLE]

where the expectation is taken over the independent samples $v_{1},\ldots,v_{N}$ from $q_{\theta_{n}}$.

2.2. Verification challenges

We now give two example model-guide pairs that illustrate verification challenges related to SVI.

The first example appears in Figure 2(a). It is the Bayesian regression example from the Pyro webpage (this example is among the benchmarks used in §8), which solves the problem of finding a line that interpolates a given set of points in $\mathbb{R}^{2}$.

The problem with this example is that the KL divergence of its model-guide pair, the main optimisation objective in SVI, is undefined. The model and guide in the figure use the random variable sigma, but they use different non-zero-probability regions, called supports, for it. In the model, the support is $[0,10]$, while that in the guide is $\mathbb{R}$. But the KL divergence from a guide to a model is defined only if for every random variable, its support in the guide is included in that in the model. We point out that this support mismatch was found by our static analyser explained in §8.

Figures 2(b) and 2(c) show two attempts to resolve the undefined-KL issue. To fix the issue, we change the distribution of sigma in the model in (b), and in the guide in (c). These revisions remove the problem about the support of sigma, but do not eliminate that of the undefined KL. In both (b) and (c), the KL divergence is $\infty$. This happens mainly because sigma can be arbitrarily close to 0 in the guide in both cases, which makes integrand in the definition of the KL divergence diverge to $\infty$.

An SVI-specific verification challenge related to this example is how to prove the well-definedness of the KL divergence and more generally the optimisation objective of an SVI algorithm. In §6.2, we provide a partial answer to the question. We give a condition for ensuring the well-definedness of the KL divergence. Our condition is more automation-friendly than the definition of KL, because it does not impose the difficult-to-check integrability requirement present in the definition of KL.

The second example appears in Figure 4(a). It uses the same model as in Figure 1(a), but has a new guide that uses a uniform distribution parameterised by $\theta\in\mathbb{R}$. For this model-guide pair, the KL divergence is well-defined for all $\theta\in\mathbb{R}$, and the optimal $\theta^{*}$ minimising the KL is $\theta^{*}=1$.

However, as shown in Figure 4(b), the gradient of the KL divergence is undefined for $\theta\in\{-1,1\}$, because the KL divergence is not differentiable at $-1$ and $1$. For all the other $\theta\in\mathbb{R}\setminus\{-1,1\}$, the KL divergence and its gradient are both defined, but the score estimator cannot estimate this gradient in an unbiased manner (i.e., in a way satisfying (1)), thereby losing the convergence guarantee to a local optimum. The precise calculation is not appropriate in this section, but we just point out that the expectation of the estimated gradient is always zero for all $\theta\in\mathbb{R}\setminus\{-1,1\}$, but the true gradient of the KL is always non-zero for those $\theta$, because it has the form: $\frac{\theta}{25}-{\mathds{1}_{[{-1\leq\theta\leq 1}]}}\cdot\frac{1}{2}\log\frac{\mathcal{N}(0;1,1)}{\mathcal{N}(0;-2,1)}.$ Here $\mathcal{N}(v;\mu,\sigma)$ is the density of the normal distribution with mean $\mu$ and standard deviation $\sigma$ (concretely, $\mathcal{N}(v;\mu,\sigma)=1/(\sqrt{2\pi}\sigma)\cdot\exp(-(v-\mu)^{2}/(2\sigma^{2}))$). The mismatch comes from the invalidity of one implicit assumption about interchanging integration and gradient in the justification of the score estimator; see §5 for detail.

To sum up, the second example shows that even if the KL divergence is defined, its gradient is sometimes undefined, and also that even if both the KL divergence and its gradient are defined, the sample-based estimate of the gradient in a standard SVI algorithm may be biased—this means that the equation similar to (1) does not hold and an SVI algorithm is no longer guaranteed to converge to a local optimum. Proving that these failure cases do not arise is another SVI-specific verification challenge. In §6.3, we give another example of similar flavour, and provide an automation-friendly condition that ensures the existence of the KL divergence and its gradient as well as the unbiasedness of the gradient estimate of the score estimator.

We conclude the section by emphasising that the aforementioned issues could have a large impact on the results of SVI. For instance, running the Bayesian regression example in Figure 2(a) in Pyro after a small change (0.4 instead of 0.05 in line 14) sometimes results in a crash. Also, running its revised version in Figure 2(c) leads to the complete failure of an optimiser and produces meaningless results. These observations strengthen the importance of resolving the verification challenges presented above.444The aforementioned issues do not always cause problems in inference results (e.g., Figures 2(a) and  2(c) mostly give reasonable results), because the random seeds and the initial values for SVI could be set well so that the probability of SVI going wrong becomes low. We emphasise, however, that the probability is non-zero.

3. Review of Measure Theory and Notations

A $\sigma$-algebra $\Sigma$ on a set $X$ is a collection of subsets of $X$ such that (i) $X\in\Sigma$; (ii) $A_{0}\cup A_{1}\in\Sigma$ and $X\setminus A_{0}\in\Sigma$ for all $A_{0},A_{1}\in\Sigma$; (iii) $\bigcup_{n}A_{n}\in\Sigma$ when all subsets $A_{n}$ are in $\Sigma$. An equivalent but easier-to-remember characterisation is that $\Sigma$ is closed under boolean operations and countable union and intersection. We call the pair $(X,\Sigma)$ of a set and a $\sigma$-algebra measurable space, and subsets in $\Sigma$ measurable. A function $f$ from a measurable space $(X,\Sigma)$ to another measurable space $(X^{\prime},\Sigma^{\prime})$ is measurable if $f^{-1}(B)\in\Sigma$ for all $B\in\Sigma^{\prime}$.

An example of measurable space is the $n$-dimensional Euclidean space $\mathbb{R}^{n}$ with the Borel $\sigma$-algebra \mathcal{B}\triangleq\sigma(\{(-\infty,r_{1})\times\cdots\times(-\infty,r_{n})~{}\big{|}~{}r\in\mathbb{R}^{n}\}), where $\sigma$ is the closure operator that converts a collection of subsets of $\mathbb{R}^{n}$ into the smallest $\sigma$-algebra containing the collection. Subsets $X$ of $\mathbb{R}^{n}$, such as $[0,\infty)^{n}$, form measurable spaces with the $\sigma$-algebra $\{X\cap A\,\mid\,A\in\mathcal{B}\}$. Another example is a set $X$ with the so called discrete $\sigma$-algebra on $X$ that consists of all subsets of $X$.

A measure $\mu$ on a measurable space $(X,\Sigma)$ is a function from $\Sigma$ to $[0,\infty]$ such that $\mu(\emptyset)=0$ and $\mu$ satisfies the countable additivity condition: for a countable family of disjoint measurable subsets $B_{n}$,

[TABLE]

A well-known example is the Lebesgue measure $\lambda^{n}$ on $\mathbb{R}^{n}$ which maps each measurable subset of $\mathbb{R}^{n}$ to its volume in the usual sense.555The Lebesgue measure $\lambda^{n}$ is the unique measure on $\mathbb{R}^{n}$ that sets the volume of the unit cube $(0,1)^{n}$ to $1$ and is translation invariant: for all measurable subsets $A$ and $r\in\mathbb{R}^{n}$, $\lambda^{n}(A)=\lambda^{n}(\{r^{\prime}-r\,\mid\,r^{\prime}\in A\})$. When $\mu(X)\leq 1$, we call $\mu$ subprobability measure. If $\mu(X)=1$, we may drop “sub”, and call $\mu$ probability measure.

The Lebesgue integral $\int$ is a partial operator that maps a measure $\mu$ on $(X,\Sigma)$ and a real-valued measurable function on the same space $(X,\Sigma)$ to a real number. It is denoted by $\int\mu(\mathrm{d}x)\,f(x)$. To follow the paper, it is enough to know that this integral generalises the usual Riemann integral from calculus.666Another useful fact is that when $f$ is non-negative, $\int\mu(\mathrm{d}x)\,f(x)=\sup\sum_{i}(\inf_{x\in A_{i}}f(x))\cdot\mu(A_{i})$ where the supremum is taken with respect to all finite partitions $\{A_{i}\}_{i}$ of $X$ into measurable subsets. For a measure $\nu$ on $(X,\Sigma)$, if $\nu(A)=\int\mu(\mathrm{d}x)\,(f(x)\cdot{\mathds{1}_{[{x\in A}]}})$ for non-negative $f$, we say that $f$ is the density of $\nu$ with respect to $\mu$ and call $\mu$ reference measure.

In the paper, we use a few well-known methods for building measurable spaces.

The first method applies when we are given a set $X$ and a collection of functions $\{f_{i}:X\to X_{i}\,\mid\,i\in I\}$ to measurable spaces $(X_{i},\Sigma_{i})$. The method is to equip $X$ with the smallest $\sigma$-algebra $\Sigma$ making all $f_{i}$’s measurable: $\Sigma\triangleq\sigma(\{f^{-1}_{i}(B)\,\mid\,i\in I,\,B\in\Sigma_{i}\})$.

The second relies on two constructions for sets, i.e., product and disjoint union. Suppose that we are given measurable spaces $(X_{i},\Sigma_{i})$ for all $i\in I$. We define a product measurable space that has $\prod_{i\in I}X_{i}$ as its underlying set and the following product $\sigma$-algebra $\bigotimes_{i\in I}\Sigma_{i}$ as its $\sigma$-algebra:

[TABLE]

The construction of the product $\sigma$-algebra can be viewed as a special case of the first method where we consider the smallest $\sigma$-algebra on $\prod_{i\in I}X_{i}$ that makes every projection map to $X_{i}$ measurable. When the $X_{i}$ are disjoint, they can be combined as disjoint union. The underlying set in this case is $\bigcup_{i\in I}X_{i}$, and the $\sigma$-algebra is

[TABLE]

When $I=\{i,j\}$ with $i\neq j$, we denote the product measurable space by $(X_{i}\times X_{j},\Sigma_{i}\otimes\Sigma_{j})$. In addition, if $X_{i}$ and $X_{j}$ are disjoint, we write $(X_{i}\cup X_{j},\Sigma_{i}\oplus\Sigma_{j})$ for the disjoint-union measurable space.

The third method builds a measurable space out of measures or a certain type of measures, such as subprobability measures. For a measurable space $(X,\Sigma)$, we form a measurable space with measures. The underlying set $\mathrm{Mea}(X)$ and $\sigma$-algebra $\Sigma_{\mathrm{Mea}(X)}$ of the space are defined by

[TABLE]

The difficult part to grasp is $\Sigma_{\mathrm{Mea}(X)}$. Once again, a good approach for understanding it is to realise that $\Sigma_{\mathrm{Mea}(X)}$ is the smallest $\sigma$-algebra that makes the function $\mu\longmapsto\mu(A)$ from $\mathrm{Mea}(X)$ to $\mathbb{R}$ measurable for all measurable subsets $A\in\Sigma$. This measurable space gives rise to a variety of measurable spaces, each having a subset $M$ of $\mathrm{Mea}(X)$ as its underlying set and the induced $\sigma$-algebra $\Sigma_{M}=\{A\cap M\,\mid\,A\in\Sigma_{\mathrm{Mea}(X)}\}$. In the paper, we use two such spaces, one induced by the set $\mathrm{Sp}(X)$ of subprobability measures on $X$ and the other by the set $\Pr(X)$ of probability measures.

A measurable function $f$ from $(X,\Sigma)$ to $(\mathrm{Mea}(Y),\Sigma_{\mathrm{Mea}(Y)})$ is called kernel. If $f(x)$ is a subprobability measure (i.e., $f(x)\in\mathrm{Sp}(Y)$) for all $x$, we say that $f$ is a subprobability kernel. In addition, if $f(x)$ is a probability measure (i.e., $f(x)\in\Pr(Y)$) for all $x$, we call $f$ probability kernel. A good heuristic is to view a probability kernel as a random function and a subprobability kernel as a random partial function. We use well-known facts that a function $f:X\to\mathrm{Mea}(Y)$ is a subprobability kernel if and only if it is a measurable map from $(X,\Sigma)$ to $(\mathrm{Sp}(Y),\Sigma_{\mathrm{Sp}(Y)})$, and that similarly a function $f$ is a probability kernel if and only if it is a measurable function from $(X,\Sigma)$ to $(\Pr(Y),\Sigma_{\Pr(Y)})$.

We use a few popular operators for constructing measures throughout the paper. We say that a measure $\mu$ on a measurable space $(X,\Sigma)$ is finite if $\mu(X)<\infty$, and $\sigma$-finite if there is a countable partition of $X$ into measurable subsets $X_{n}$’s such that $\mu(X_{n})<\infty$ for every $n$. Given a finite or countable family of $\sigma$-finite measures $\{\mu_{i}\}_{i\in I}$ on measurable spaces $(X_{i},\Sigma_{i})$’s, the product measure of $\mu_{i}$’s, denoted $\bigotimes_{i\in I}\mu_{i}$, is the unique measure on $(\prod_{i\in I}X_{i},\bigotimes_{i\in I}\Sigma_{i})$ such that for all measurable subsets $A_{i}$ of $X_{i}$,

[TABLE]

Given a finite or countable family of measures $\{\mu_{i}\}_{i\in I}$ on disjoint measurable spaces $(X_{i},\Sigma_{i})$’s, the sum measure of $\mu_{i}$’s, denoted $\bigoplus_{i\in I}\mu_{i}$, is the unique measure on $(\sum_{i\in I}X_{i},\bigoplus_{i\in I}\Sigma_{i})$ such that

[TABLE]

Throughout the paper, we take the convention that the set $\mathbb{N}$ of natural numbers includes [math]. For all positive integers $n$, we write $[n]$ to mean the set $\{1,2,\ldots,n\}$.

4. Simple Probabilistic Programming Language

In this section, we describe the syntax and semantics of a simple probabilistic programming language, which we use to present the theoretical results of the paper. The measure semantics in §4.2 uses results and observations from (Staton et al., 2016), but the density semantics and the other materials in §4.3 are new in this work.

4.1. Syntax

We use an extension of the standard while language with primitives for probabilistic programming. The grammar of the language is given in Figure 5. Variables in the language store real numbers, but expressions may denote reals, booleans and strings, and they are classified into $E,B,S$ based on these denoted values. The primitive functions $f$ for reals and $g$ for strings may be usual arithmetic and string operations, such as multiplication and exponentiation for $f$ and string concatenation for $g$.

The grammar for $C$ includes the cases for the standard constructs of the while language, such as assignment, sequencing, conditional statement and loops. In addition, it has two constructs for probabilistic programming. The first $x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})$ draws a sample from the normal distribution with mean $E_{1}$ and standard deviation $E_{2}$ and names the sample with the string $S$. The next $\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})$ expresses that a sample is drawn from the normal distribution with mean $E_{1}$ and standard deviation $E_{2}$ and the value of this sample is observed to be $E_{0}$. It lets the programmers express information about observed data inside programs. Operationally, this construct can be understood as an instruction for updating a global variable that stores the so called importance score of the execution. The score quantitatively records how well the random choices in the current execution match the observations, and the score statement updates this score by multiplying it with the density at $E_{0}$ of the appropriate normal distribution.

Consider the following program:

[TABLE]

The program specifies a model with one random variable $``a"$. Using a relatively flat normal distribution, the program specifies a prior belief that the value of the random variable $``a"$ is likely to be close to $0.0$ and lie between $-2\times 5.0$ and $2\times 5.0$. The next score statement refines this belief with one data point $3.0$, which is a noisy observation of the value of $``a"$ (bound to $x$). The parameters to the normal density in the statement express that the noise is relatively small, between $-2\times 1.0$ and $2\times 1.0$. Getting the refined belief, called posterior distribution, is the reason that a data scientist writes a model like this program. It is done by an inference algorithm of the language.

Permitting only the normal distribution does not limit the type of models expressible in the language. Every distribution can be obtained by transforming the standard normal distribution.777Here we consider only Borel spaces. Although using only the normal distribution does not affect the expressiveness, it has an impact on stochastic variational inference to be discussed later, because it requires a guide to use only normal distributions. But the impact is not significant, because most well-known approaches for creating guides from the machine-learning literature (such as extensions of variational autoencoder) use normal distributions only, or can be made to do so easily.

4.2. Measure semantics

The denotational semantics of the language just presented is mostly standard, but employs some twists to address the features for probabilistic programming (Staton et al., 2016).

Here is a short high-level overview of the measure semantics. Our semantics defines multiple measurable spaces, such as $\mathit{{Store}}$ and $\mathit{{State}}$, that hold mathematical counterparts to the usual actors in computation, such as program stores (i.e., mappings from variables to values) and states (which consist of a store and further components). Then, the semantics interprets expressions $E,B,S$ and commands $C$ as measurable functions of the following types:

[TABLE]

Here $\mathbb{R}$ is the measurable space of reals with the Borel $\sigma$-algebra, and $\mathbb{B}$ and $\mathit{{Str}}$ are discrete measurable spaces of booleans and strings. $\mathit{{Store}}$ and $\mathit{{State}}$ are measurable spaces for stores (i.e., maps from variables to values) and states which consist of a store and a part for recording information about sampled random variables. Note that the target measurable space of commands is built by first taking the product of measurable spaces $\mathit{{State}}$ and $[0,\infty)$ and then forming a space out of subprobability measures on $\mathit{{State}}\times[0,\infty)$. This construction indicates that commands denote probabilistic computations, and the result of each such computation consists of an output state and a score which expresses how well the computation matches observations expressed with the score statements in the command $C$. Some of the possible outcomes of the computation may lead to non-termination or an error, and these abnormal outcomes are not accounted for by the semantics, which is why ${\llbracket C\rrbracket}(\sigma)$ for a state $\sigma\in\Sigma$ is a subprobability distribution. The semantics of expressions is much simpler. It just says that expressions do not involve any probabilistic computations, so that they denote deterministic measurable functions.

We now explain how this high-level idea gets implemented in our semantics. Let $\mathit{{Var}}$ be a countably infinite set of variables. Our semantics uses the following sets:

[TABLE]

where $[X\to Y]$ denotes the set of all functions from $X$ to $Y$. A state $\sigma$ consists of a store $s$ and a random database $r$. The former fixes the values of variables, and the latter records the name (given as a string) and the value of each sampled random variable. The domain of $r$ is the names of all the sampled random variables. By insisting that $r$ should be a map, the semantics asserts that no two random variables have the same name. For each state $\sigma$, we write $\sigma_{s}$ and $\sigma_{r}$ for its store and random database components, respectively. Also, for a variable $x$, a string $\alpha$ and a value $v$, we write $\sigma[x\mapsto v]$ and $\sigma[\alpha\mapsto v]$ to mean $(\sigma_{s}[x\mapsto v],\sigma_{r})$ and $(\sigma_{s},\sigma_{r}[\alpha\mapsto v])$.

We equip all of these sets with $\sigma$-algebras and turn them to measurable spaces in a standard way. Note that we constructed the sets from $\mathbb{R}$ by repeatedly applying the product and disjoint-union operators. We equip $\mathbb{R}$ with the usual Borel $\sigma$-algebra. Then, we parallel each usage of the product and the disjoint-union operators on sets with that of the corresponding operators on $\sigma$-algebras. This gives the $\sigma$-algebras for all the sets defined above. Although absent from the above definition, the measurable spaces $\mathbb{B}$ and $\mathit{{Str}}$ equipped with discrete $\sigma$-algebras are also used in our semantics.

We interpret expressions $E,B,S$ as measurable functions ${\llbracket E\rrbracket}:\mathit{{Store}}\to\mathbb{R}$, ${\llbracket B\rrbracket}:\mathit{{Store}}\to\mathbb{B}$, and ${\llbracket S\rrbracket}:\mathit{{Store}}\to\mathit{{Str}}$, under the assumption that the semantics of primitive real-valued $f$ of arity $n$ and string-valued $g$ of arity $(m,l)$ are given by measurable functions ${\llbracket f\rrbracket}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ and ${\llbracket g\rrbracket}:\mathit{{Str}}^{m}\times\mathbb{R}^{l}\rightarrow\mathit{{Str}}$. It is standard, and we describe it only for some sample cases of $E$ and $S$:

[TABLE]

Lemma 4.1.

For all expressions $E$, $B$, and $S$, their semantics ${\llbracket E\rrbracket}$, ${\llbracket B\rrbracket}$ and ${\llbracket S\rrbracket}$ are measurable functions from $\mathit{{Store}}$ to $\mathbb{R}$, $\mathbb{B}$ and $\mathit{{Str}}$, respectively.

We interpret commands $C$ as measurable functions from $\mathit{{State}}$ to $\mathrm{Sp}(\mathit{{State}}\times[0,\infty))$, i.e., subprobability kernels from $\mathit{{State}}$ to $\mathit{{State}}\times[0,\infty)$. Let $\mathcal{K}$ be the set of subprobability kernels from $\mathit{{State}}$ to $\mathit{{State}}\times[0,\infty)$, and $\Sigma_{\mathit{{State}}\times[0,\infty)}$ be the $\sigma$-algebra of the product space ${\mathit{{State}}\times[0,\infty)}$. We equip $\mathcal{K}$ with the partial order $\sqsubseteq$: for all $\kappa,\kappa^{\prime}\in\mathcal{K}$, $\kappa\sqsubseteq\kappa^{\prime}$ if and only if $\kappa(\sigma)(A)\leq\kappa^{\prime}(\sigma)(A)$ for all $\sigma\in\mathit{{State}}$ and $A\in\Sigma_{\mathit{{State}}\times[0,\infty)}$. The next lemma is a minor adaptation of a known result (Kozen, 1981).

Lemma 4.2.

$(\mathcal{K},\sqsubseteq)$* is an $\omega$-complete partial order with the least element $\bot\triangleq(\lambda\sigma.\,\lambda A.\,0)$.*

The measure semantics of a command $C$ is defined in Figure 6.

The interpretation of the loop is the least fixed point of the function $F$ on the $\omega$-complete partial order $\mathcal{K}$. The function $F$ is continuous, so the least fixed point is obtained by the $\omega$-limit of the sequence $\{F^{n}(\bot)\}_{n}$. In the definition of $F$, the argument $\kappa$ is a sub-probability kernel, and it represents the computation after the first iteration of the loop. The semantics of the sample statement uses an indicator function to exclude erroneous executions where the argument $S$ denotes a name already used by some previous random variable, or the standard deviation ${\llbracket E_{2}\rrbracket}\sigma_{s}$ is not positive. When this check passes, it distils $A$ to a property on the value of $x$ and computes the probability of the property using the normal distribution with mean ${\llbracket E_{1}\rrbracket}\sigma_{s}$ and standard deviation ${\llbracket E_{2}\rrbracket}\sigma_{s}$.

Theorem 4.3.

For every command $C$, its interpretation ${\llbracket C\rrbracket}$ is well-defined and belongs to $\mathcal{K}$.

4.3. Posterior inference and density semantics

We write a probabilistic program to answer queries about the model and data that it describes. Among such queries, posterior inference is one of the most important and popular. Let ${\sigma_{I}}=(s_{I},[])$ be the initial state that consists of some fixed store and the empty random database. In our setting, posterior inference amounts to finding information about the following probability measure $\Pr(C,\cdot)$ for a command $C$. For a measurable $A\subseteq\mathit{{RDB}}$,

[TABLE]

The probability measure $\Pr(C,\cdot)$ is called the posterior distribution of $C$, and $\mathrm{Mea}(C,\cdot)$ the unnormalised posterior distribution of $C$ (in $\mathrm{Mea}(C,\cdot)$ and $\Pr(C,\cdot)$, we elide the dependency on $s_{I}$ to avoid clutter). Finding information about the former is the goal of most inference engines of existing probabilistic programming languages. Of course, $\Pr(C,\cdot)$ is not defined when the normalising constant $Z_{C}$ is infinite or zero. The inference engines regard such a case as an error that a programmer should avoid, and consider only $C$ without those errors.

Most algorithms for posterior inference use the density semantics of commands. They implicitly pick measures on some measurable spaces used in the semantics. These measures are called reference measures, and constructed out of Lebesgue and counting measures (Bhat et al., 2012, 2013; Hur et al., 2015). Then, the algorithms interpret commands as density functions with respect to these measures. One outcome of this density semantics is that the unnormalised posterior distribution $\mathrm{Mea}(C,\cdot)$ of a command $C$ has a measurable function $f:\mathit{{RDB}}\to[0,\infty)$ such that $\mathrm{Mea}(C,A)=\int\rho(\mathrm{d}r)\,({\mathds{1}_{[{r\in A}]}}\cdot f(r))$, where $\rho$ is a reference measure on $\mathit{{RDB}}$. Function $f$ is called density of $\mathrm{Mea}(C,\cdot)$ with respect to $\rho$.

In the rest of this subsection, we reformulate the semantics of commands using density functions. To do this, we need to set up some preliminary definitions.

First, we look at a predicate and an operator for random databases, which are about the possibility and the very act of merging two databases. For $r,r^{\prime}\in\mathit{{RDB}}$, define the predicate $r\#r^{\prime}$ by:

[TABLE]

When $r\#r^{\prime}$, let $r\uplus r^{\prime}$ be the random database obtained by merging $r$ and $r^{\prime}$:

[TABLE]

Lemma 4.4.

For every measurable $h:\mathit{{RDB}}\times\mathit{{RDB}}\times\mathit{{RDB}}\to\mathbb{R}$, the function $(r,r^{\prime})\longmapsto{\mathds{1}_{[{r\#r^{\prime}}]}}\times h(r,r^{\prime},r\uplus r^{\prime})$ from $\mathit{{RDB}}\times\mathit{{RDB}}$ to $\mathbb{R}$ is measurable.

Second, we define a reference measure $\rho$ on $\mathit{{RDB}}$:

[TABLE]

where $\lambda$ is the Lebesgue measure on $\mathbb{R}$. As explained in the preliminary section, the symbol $\otimes$ here represents the operator for constructing a product measure. In particular, $\bigotimes_{\alpha\in K}\lambda$ refers to the product of the $|K|$ copies of the Lebesgue measure $\lambda$ on $\mathbb{R}$. In the above definition, we view functions in $[K\to\mathbb{R}]$ as tuples with $|K|$ real components and measure sets of such functions using the product measure $\bigotimes_{\alpha\in K}\lambda$. When $K$ is the empty set, $\bigotimes_{\alpha\in K}\lambda$ is the nullary-product measure on $\{[]\}$, which assigns $1$ to $\{[]\}$ and [math] to the empty set.

The measure $\rho$ computes the size of each measurable subset $R\subseteq\mathit{{RDB}}$ in three steps. It splits a given $R$ into groups based on the domains of elements in $R$. Then, it computes the size of each group separately, using the product of the Lebesgue measure. Finally, it adds the computed sizes. The measure $\rho$ is not finite, but it satisfies the $\sigma$-finiteness condition,888The condition lets us use Fubini theorem when showing the well-formedness of the density semantics in this subsection and relating this semantics with the measure semantics in the previous subsection. the next best property.

Third, we define a partially-ordered set $\mathcal{D}$ with certain measurable functions. We say that a function $g:\mathit{{Store}}\times\mathit{{RDB}}\to\{\bot\}\cup(\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty))$ uses random databases locally or is local if for all $s,s^{\prime}\in\mathit{{Store}}$, $r,r^{\prime}\in\mathit{{RDB}}$, and $w^{\prime},p^{\prime}\in[0,\infty)$,

[TABLE]

The condition describes the way that $g$ uses a given random database $r$, which plays the role of a bank of random seeds (that $g$ may partially consume as it needs random values). Some part $r^{\prime\prime}$ of $r$ may be consumed by $g$, but the unconsumed part $r^{\prime}$ of $r$ does not change and is returned in the output. Also, the behaviour of $g$ does not depend on the unconsumed $r^{\prime}$. We define the set $\mathcal{D}$ by

[TABLE]

Here we view $\{\bot\}$ and $[0,\infty)$ as measurable spaces equipped with discrete and Borel $\sigma$-algebras. Also, we regard $\mathit{{Store}}\times\mathit{{RDB}}$ and $\{\bot\}\cup(\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty))$ as measurable spaces constructed by the product and disjoint-union operators on measurable spaces, as explained in §3.

The locality in the definition of $\mathcal{D}$ formalises expected behaviours of commands. In fact, as we will show shortly, it is satisfied by all commands in our density semantics. This property plays an important role when we establish the connection between the density semantics in this subsection and the standard measure semantics in §4.2.

The functions in $\mathcal{D}$ are ordered pointwise: for all $g,g^{\prime}\in\mathcal{D}$,

[TABLE]

Lemma 4.5.

$(\mathcal{D},\sqsubseteq)$* is an $\omega$-complete partial order and has the least element $a\longmapsto\bot$. Thus, every continuous function $G$ on $\mathcal{D}$ has a least fixed point (and this least fixed point is unique).*

For each $g\in\mathcal{D}$, let $g^{\ddagger}$ be the following lifting to a function on $\{\bot\}\cup(\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty))$:

[TABLE]

This lifting lets us compose two functions in $\mathcal{D}$.

Using these preliminary definitions, we define a density semantics in Figure 7, where a command $C$ means a function ${\llbracket C\rrbracket}_{d}\in\mathcal{D}$. The notation $r\setminus v$ in the figure means the removal of the entry $v$ from the finite map $r$ if $v\in\mathrm{dom}(r)$; otherwise, $r\setminus v$ is just $r$. The set membership ${\llbracket C\rrbracket}_{d}\in\mathcal{D}$ says that ${\llbracket C\rrbracket}_{d}$ is a local measurable function from $\mathit{{Store}}\times\mathit{{RDB}}$ to $\{\bot\}\cup\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty)$. Thus, the function ${\llbracket C\rrbracket}_{d}$ takes a store $s$ and a random database $r$ as inputs, where the former fixes the values of variables at the start of $C$ and the latter specifies random seeds some of which $C$ may consume to sample random variables. Given such inputs, the function outputs an updated store $s^{\prime}$, the part $r^{\prime}$ of $r$ not consumed by $C$, the total score $w^{\prime}$ expressing how well the execution of $C$ matches observations, and the probability density $p^{\prime}$ of $C$ at the consumed part of $r$. If $r$ does not contain enough random seeds, or the execution of $C$ encounters some runtime error, or it falls into an infinite loop, then the function returns $\bot$.

Lemma 4.6.

For every command $C$, its semantics ${\llbracket C\rrbracket}_{d}$ is well-defined and belongs to $\mathcal{D}$.

The density semantics ${\llbracket C\rrbracket}_{d}$ is closely related to the measure semantics ${\llbracket C\rrbracket}$ defined in §4.2. Both interpretations of $C$ describe the computation of $C$ but from slightly different perspectives. To state this relationship formally, we need a few notations. For $g\in\mathcal{D}$ and $s\in\mathit{{Store}}$, define

[TABLE]

Both $\mathit{{dens}}(g,s)$ and $\mathit{{get}}(g,s)$ are concerned with random databases $r$ that precisely describe the randomness needed by the execution of $g$ from $s$. This is formalised by the use of $[]$ in the definitions. The function $\mathit{{dens}}(g,s)$ assigns a score to such an $r$, and in so doing, it defines a probability density on $\mathit{{RDB}}$ with respect to the reference measure $\rho$. The function $\mathit{{get}}(g,s)$ computes a store that the computation of $g$ would result in when started with such an $r$. We often write $\mathit{{dens}}(C,s)$ and $\mathit{{get}}(C,s)$ to mean $\mathit{{dens}}({\llbracket C\rrbracket}_{d},s)$ and $\mathit{{get}}({\llbracket C\rrbracket}_{d},s)$, respectively.

Lemma 4.7.

For all $g\in\mathcal{D}$, the following functions from $\mathit{{Store}}\times\mathit{{RDB}}$ to $\mathbb{R}$ and $\{\bot\}\cup\mathit{{Store}}$ are measurable: $(s,r)\longmapsto\mathit{{dens}}(g,s)(r)$ and $(s,r)\longmapsto\mathit{{get}}(g,s)(r)$.

The next lemma is the main reason that we considered the locality property. It plays a crucial role in proving Theorem 4.9, the key result of this subsection.

Lemma 4.8.

For all non-negative bounded measurable functions $h:(\{\bot\}\cup\mathit{{Store}})\times\mathit{{RDB}}\to\mathbb{R}$, stores $s$, and functions $g_{1},g_{2}\in\mathcal{D}$, we have that

[TABLE]

Assume that $(\{\bot\}\cup\mathit{{Store}})\times\mathit{{RDB}}$ is ordered as follows: for all $(a,r),(a^{\prime},r^{\prime})\in(\{\bot\}\cup\mathit{{Store}})\times\mathit{{RDB}}$,

[TABLE]

Theorem 4.9.

For all non-negative bounded measurable monotone functions $h:(\{\bot\}\cup\mathit{{Store}})\times\mathit{{RDB}}\to\mathbb{R}$ and states $\sigma$,

[TABLE]

Corollary 4.10.

$\mathrm{Mea}(C,A)=\int\rho(\mathrm{d}r)\,({\mathds{1}_{[{r\in A}]}}\cdot\mathit{{dens}}(C,s_{I})(r))$* for all $C$ and all measurable $A$.*

Proof.

We instantiate Theorem 4.9 with $h(a,r)={\mathds{1}_{[{r\in A}]}}$ and $\sigma=\sigma_{I}$. Recall that $(\sigma_{I})_{s}=s_{I}$ and $(\sigma_{I})_{r}=[]$. Thus, the conclusion of Theorem 4.9 in this case says that

[TABLE]

This gives the equation claimed by the corollary. ∎

This corollary says that $\mathit{{dens}}(C,s_{I})$ is the density of the measure $\mathrm{Mea}(C,\cdot)$ with respect to $\rho$, and supports our informal claim that ${\llbracket C\rrbracket}_{d}$ computes the density of the measure ${\llbracket C\rrbracket}$.

5. Stochastic variational inference

In this section, we explain stochastic variational inference (SVI) algorithms using the semantics that we have developed so far. In particular, we describe the requirements made implicitly by one fundamental SVI algorithm, which is regarded most permissive by the ML researchers because the algorithm does not require the differentiability of the density of a given probabilistic model.

We call a command $C$ model if it has a finite nonzero normalising constant:

[TABLE]

where ${\sigma_{I}}$ is the initial state. Given a model $C$, the SVI algorithms attempt to infer a good approximation of $C$’s posterior distribution $\Pr(C,\cdot)$ defined in (2). They tackle this posterior-inference problem in two steps.

First, the SVI algorithms fix a collection of approximate distributions. They usually do so by asking the developer of $C$ to provide a command $D_{\theta}$ parameterised by $\theta\in\mathbb{R}^{p}$, which can serve as a template for approximation distributions. The command $D_{\theta}$ typically has a control-flow structure similar to that of $C$, but it is simpler than $C$: it does not use any score statements, and may replace complex computation steps of $C$ by simpler ones. In fact, $D_{\theta}$ should satisfy two formal requirements, which enforce this simplicity. The first is

[TABLE]

which means that the normalising constant of $D_{\theta}$ is $1$. The second is that $D_{\theta}$ should keep the score (i.e., the $w$ component) to be $1$, i.e.,

[TABLE]

Meeting these requirements is often not too difficult. A common technique is to ensure that $D_{\theta}$ does not use the score statement and always terminates. Figure 8 gives an example of $(C,D_{\theta})$ for simple Bayesian linear regression with three data points. Note that in this case, $D_{\theta}$ is obtained from $C$ by deleting the score statements and replacing the arguments $0.0$ and $5.0$ of normal distributions by parameter $\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})$. Following the terminology of Pyro, we call a parameterised command $D_{\theta}$ guide if it satisfies the two requirements just mentioned.

Second, the SVI algorithms search for a good parameter $\theta$ that makes the distribution described by $D_{\theta}$ close to the posterior of $C$. Concretely, they formulate an optimisation problem where the optimisation objective expresses that some form of distance from $D_{\theta}$’s distribution to $C$’s posterior should be minimised. Then, they solve the problem by a version of gradient descent.

The KL divergence is a standard choice for distance. Let $\mu,\mu^{\prime}$ be measures on $\mathit{{RDB}}$ that have densities $g$ and $g^{\prime}$ with respect to the measure $\rho$. The KL divergence from $g$ to $g^{\prime}$ is defined by

[TABLE]

In words, it is the ratio of densities $g$ and $g^{\prime}$ averaged according to $g$. If $g=g^{\prime}$, the ratio is always $1$, so that the KL becomes [math]. The KL divergence is defined only if the following conditions are met:

• •

Absolute continuity: $g^{\prime}(r)=0\implies g(r)=0$ for all $r\in\mathit{{RDB}}$,999This condition can be relaxed in a more general formulation of the KL divergence stated in terms of the so called Radon-Nikodym derivative. We do not use the relaxed condition to reduce the amount of materials on measure theory in the paper. which ensures that the integrand in (4) is well-defined even when the denominator $g^{\prime}(r)$ in (4) takes the value [math];

• •

Integrability: the integral in (4) has a finite value.

Using our semantics, we can express the KL objective as follows:

[TABLE]

Recall that $\mathit{{dens}}(D_{\theta},s_{I})$ and $\mathit{{dens}}(C,s_{I})/Z_{C}$ are densities of the probability measures of the command $D_{\theta}$ and the posterior of $C$ (Corollary 4.10), and they are defined by means of our density semantics in §4.3. Most SVI engines solve this optimisation problem by a version of gradient descent.

In the paper, we consider one of the most fundamental and versatile SVI algorithms. The algorithm is called score estimator or REINFORCE, and it works by estimating the gradient of the objective in (5) using samples and performing the gradient descent with this estimated gradient. More concretely, the algorithm starts by initialising $\theta$ with some value (usually chosen randomly) and updating it repeatedly by the following procedure:

[TABLE]

Here $N$ and $\eta$ are hyperparameters to this algorithm, the former determining the number of samples used to estimate the gradient and the latter, called learning rate, deciding how much the algorithm should follow the direction of the estimated gradient. Although we do not explain here, sampling $r_{1},\ldots,r_{N}$ and computing all of $\mathit{{dens}}(D_{\theta},s_{I})(r_{i})$, $\mathit{{dens}}(C,s_{I})(r_{i})$ and $\nabla_{\theta}(\log\mathit{{dens}}(D_{\theta},s_{I})(r_{i}))$ can be done by executing $D_{\theta}$ and $C$ multiple times under slightly unusual operational semantics (Yang, 2019). The SVI engines of Pyro and Anglican implement such operational semantics.

The average over the $N$ terms in the $\theta$-update step is the core of the algorithm. It approximates the gradient of the optimisation objective in (5):

[TABLE]

The average satisfies an important property called unbiasedness, summarised by Theorem 5.1.

Theorem 5.1.

*Let $C$ be a model, $D_{\theta}$ be a guide, and $N\neq 0\in\mathbb{N}$. Define $\mathrm{KL}_{(-)}:\mathbb{R}^{p}\to\mathbb{R}_{\geq 0}$ as $\mathrm{KL}_{\theta}$ $\triangleq$ $\mathrm{KL}(\mathit{{dens}}(D_{\theta},s_{I}){||}\mathit{{dens}}(C,s_{I})/Z_{C})$. Then, $\mathrm{KL}_{(-)}$ is well-defined and continuously differentiable with *

[TABLE]

if

1. (R1)

$\mathit{{dens}}(C,s_{I})(r)=0\implies\mathit{{dens}}(D_{\theta},s_{I})(r)=0$, for all $r\in\mathit{{RDB}}$ and $\theta\in\mathbb{R}^{p}$; 2. (R2)

for all $(r,\theta,j)\in\mathit{{RDB}}\times\mathbb{R}^{p}\times[p]$, the function $v\longmapsto\mathit{{dens}}(D_{\theta[j:v]},s_{I})(r)$ on $\mathbb{R}$ is differentiable; 3. (R3)

for all $\theta\in\mathbb{R}^{p}$,

[TABLE] 4. (R4)

for all $(\theta,j)\in\mathbb{R}^{p}\times[p]$, the function

[TABLE]

on $\mathbb{R}$ is continuously differentiable; 5. (R5)

for all $\theta\in\mathbb{R}^{p}$,

[TABLE] 6. (R6)

for all $\theta\in\mathbb{R}^{p}$,

[TABLE]

Here $\theta[j:v]$ denotes a vector in $\mathbb{R}^{p}$ that is the same as $\theta$ except that its $j$-th component is $v$.

The conclusion of this theorem (Equation (6)) and its proof are well-known (Ranganath et al., 2014), but the requirements in the theorem (and the continuous differentiability of $\mathrm{KL}_{\theta}$ in the conclusion) are rarely stated explicitly in the literature.

The correctness of the algorithm crucially relies on the unbiasedness property in Theorem 5.1. The property ensures that the algorithm converges to a local minimum with probability $1$. Thus, it is important that the requirements in the theorem are met. In fact, some of the requirements there are needed even to state the optimisation objective in (5), because without them, the objective does not exist. In the next section, we describe conditions that imply those requirements and can serve as target properties of program analysis for probabilistic programs. The latter point is worked out in detail in §7 and §8 where we discuss program analysis for probabilistic programs and SVI.

6. Conditions for stochastic variational inference

Ideally we want to have program analysers that discharge the six requirements R1-R6 in Theorem 5.1. However, except R1 and R2, the requirements are not ready for serving as the targets of static analysis algorithms. Automatically discharging them based on the first principles (such as the definition of integrability with respect to a measure) may be possible, but seems less immediate than doing so using powerful theorems from continuous mathematics.

In this section, we explain conditions that imply the requirements R3-R6 and are more friendly to program analysis than the requirements themselves. The conditions are given in two boxes (13) and (17). Throughout the section, we fix a model $C$ and a guide $D_{\theta}$.

6.1. Assumption

Throughout the section, we assume that the densities of $D_{\theta}$ and $C$ have the following form:

[TABLE]

where

• •

$M,N_{i}\in\mathbb{N}\setminus\{0\}$;

• •

$A_{1},\ldots,A_{M}$ are disjoint measurable subsets of $\mathit{{RDB}}$;

• •

$K_{i}$’s are finite sets of strings such that $\mathrm{dom}(r)=K_{i}$ for all $r\in A_{i}$;

• •

$\mu_{(i,\alpha)}$ and $\sigma_{(i,\alpha)}$ are measurable functions from $\mathbb{R}^{p}$ to $\mathbb{R}$ and $(0,\infty)$, respectively;

• •

$\mu^{\prime}_{(i,\alpha)}$ and $\mu^{\prime\prime}_{(i,j)}$ are measurable functions from $[K_{i}\to\mathbb{R}]$ to $\mathbb{R}$;

• •

$\sigma^{\prime}_{(i,\alpha)}$ and $\sigma^{\prime\prime}_{(i,j)}$ are measurable functions from $[K_{i}\to\mathbb{R}]$ to $(0,\infty)$;

• •

$c_{(i,j)}$ is a real number.

In programming terms, our assumption first implies that both $C$ and $D_{\theta}$ use at most a fixed number of random variables. That is, the number of random variables they generate must be finite not only within a single execution but also across all possible executions since the names of all random variables are found in a finite set $\bigcup_{i=1}^{M}K_{i}$. This property is met if the number of steps in every execution of $C$ and $D_{\theta}$ from $\sigma_{I}$ is bounded by some $T$ and the executions over the same program path use the same set of random variables. Note that the bound may depend on $s_{I}$. Such a bound exists for most probabilistic programs.101010Notable exceptions are models using probabilistic grammars or those from Bayesian nonparametrics. Also, the assumption says that the parameters of normal distributions in sample statements in $D_{\theta}$ may depend only on $\theta$, but not on other sampled random variables. This is closely related to a common approach for designing approximate distributions in variational inference, called mean-field approximation, where the approximate distribution consists of independent normal random variables.

We use the term “assumption” here, instead of “condition” in the following subsections because the assumed properties are rather conventional and they are not directly related to the requirements in Theorem 5.1, at least not as much as the conditions that we will describe next.

6.2. Condition for the requirement R3

Note that the integral in the requirement R3 can be written as the sum of two expectations:

[TABLE]

The minus of the first term (i.e., $-\mathbb{E}_{\mathit{{dens}}(D_{\theta},s_{I})(r)}\left[{\log\mathit{{dens}}(D_{\theta},s_{I})(r)}\right]$) is called the differential entropy of the density $\mathit{{dens}}(D_{\theta},s_{I})$. Intuitively, it is large when the density on $\mathbb{R}^{n}$ is close to the Lebesgue measure, which is regarded to represent the absence of information. The differential entropy is sometimes undefined (Ghourchian et al., 2017). Fortunately, a large class of probability densities (containing many commonly used probability distributions) have well-defined entropies (Nair et al., 2006; Ghourchian et al., 2017). Our $\mathit{{dens}}(D_{\theta},s_{I})$ is one of such fortunate cases.

Theorem 6.1.

$\mathbb{E}_{\mathit{{dens}}(D_{\theta},s_{I})(r)}\left[{|\log\mathit{{dens}}(D_{\theta},s_{I})(r)|}\right]<\infty$* under our assumption in §6.1.*

We remark that a violation of the assumption (7) for $D_{\theta}$ can result in an undefined entropy, as illustrated by the following examples.

Example 6.2.

Consider guides $D_{(i,\theta)}$ defined as follows ($i=1,2$):

[TABLE]

where for some $n\geq 1$ and $c\neq 0\in\mathbb{R}$,111111Formally, we should implement $E_{1}[x_{1}]$ as an application of a primitive function $f_{1}$ to $x_{1}$ that has the semantics described by the if-then-else statement here.

[TABLE]

None of $\mathit{{dens}}(D_{(i,\theta)},s_{I})$’s satisfies the assumption (7) because the standard deviation of the normal distribution for $x_{2}$ depends on the value of $x_{1}$. The entropies of $\mathit{{dens}}(D_{(i,\theta)},s_{I})$’s are all undefined: $\mathbb{E}_{\mathit{{dens}}(D_{(i,\theta)},s_{I})(r)}{[|\log\mathit{{dens}}(D_{(i,\theta)},s_{I})(r)|]}=\infty$ for all $i=1,2$. ∎

Since the first term of (8) is always finite by Theorem 6.1, it is enough to ensure the finiteness of the second term of (8). For $i\in[M]$, define the set of (absolute) affine functions on $[K_{i}\to\mathbb{R}]$ as:

[TABLE]

Our condition for ensuring the finiteness of the second term is as follows:

[TABLE]

Theorem 6.3.

The condition (13) implies $\mathbb{E}_{\mathit{{dens}}(D_{\theta},s_{I})(r)}\left[{|\log\mathit{{dens}}(C,s_{I})(r)|}\right]<\infty$ under our assumption in §6.1. Thus, in that case, it entails the requirement R3 (i.e., the objective in (5) is well-defined).

Our condition in (13) is sufficient but not necessary for the objective in (5) to be defined. However, its violation is a good warning, as illustrated by our next examples.

Example 6.4.

Consider models $C_{1},\ldots,C_{4}$ and a guide $D_{\theta}$ defined as follows:

[TABLE]

where for some $n\geq 1$ and $c\neq 0\in\mathbb{R}$,

[TABLE]

Let $A=[\{``a_{1}",``a_{2}"\}\to\mathbb{R}]$. For $r\in A$, define

[TABLE]

Then, we have that

[TABLE]

None of $\mu^{\prime}_{1}$, $\mu^{\prime}_{2}$, $\sigma^{\prime}_{3}$, and $\sigma^{\prime}_{4}$ satisfies the condition in (13). The function $\mu^{\prime}_{1}$ is not bounded in $\{r\in A\,\mid\,-1\leq r(``a_{1}")\leq 1\wedge-1\leq r(``a_{2}")\leq 1\}$, but every $\mu^{\prime}$ satisfying the condition in (13) should be bounded. Also, the cubic exponential growth of $\mu^{\prime}_{2}$ cannot be bounded by any linear exponential function on $r$. The violation of the condition by $\sigma^{\prime}_{3}$ and $\sigma^{\prime}_{4}$ can be shown similarly.

In fact, the objective in (5) is not defined for all of the four cases. This is because for all $i=1,\ldots,4$, $\mathbb{E}_{\mathit{{dens}}(D_{\theta},s_{I})(r)}\left[{|\log\mathit{{dens}}(C_{i},s_{I})(r)|}\right]=\infty$. ∎

We now show that the condition in (13) is satisfied by a large class of functions in machine learning applications, including functions using neural networks. We call a function $\mathit{nn}:\mathbb{R}^{n}\to\mathbb{R}$ an affine-bounded neural network if there exist functions $f_{j}:\mathbb{R}^{n_{j}}\to\mathbb{R}^{n_{j+1}}$ and affine functions $l_{j}:\mathbb{R}^{n_{j}}\to\mathbb{R}$ for all $1\leq j\leq d$ such that (i) $n_{1}=n$ and $n_{d+1}=1$; (ii) $\mathit{{nn}}=f_{d}\circ\cdots\circ f_{1}$; and (iii) $\|f_{j}(v)\|_{1}\leq l_{j}(|v_{1}|,\ldots,|v_{n_{j}}|)$ for all $1\leq j\leq d$ and $v\in\mathbb{R}^{n_{j}}$, where $\|\cdot\|_{1}$ denotes the $\ell_{1}$-norm. Note that each component of $f_{j}$ can be, for instance, an affine function, one of commonly used activation functions (e.g., relu, tanh, sigmoid, softplus), or one of min/max functions ($\min$ and $\max$). Therefore, most of neural networks used in machine learning applications are indeed affine-bounded. Lemma 6.5 indicates that a wide range of functions satisfy the condition (13).

Lemma 6.5.

Pick $i\in[M]$ and $\alpha\in K_{i}$. Let $(\alpha_{1},\ldots,\alpha_{J})$ be an enumeration of the elements in $K_{i}$, and $\overline{r}\triangleq(r(\alpha_{1}),\ldots,r(\alpha_{J}))\in\mathbb{R}^{J}$ be an enumeration of the values of $r\in A_{i}$. Consider any affine-bounded neural network $\mathit{nn}:\mathbb{R}^{J}\to\mathbb{R}$, polynomial $\mathit{poly}:\mathbb{R}^{J}\to\mathbb{R}$, and $c\in(0,\infty)$. Then, the below list of functions $\mu^{\prime}_{(i,\alpha)}$ and $\sigma^{\prime}_{(i,\alpha)}$ on $A_{i}$ satisfy the condition in (13):

[TABLE]

where $\mathrm{softplus}(v)\triangleq\log(1+\exp(v))$. Moreover, the same holds for $\mu^{\prime\prime}_{(i,j)}$ and $\sigma^{\prime\prime}_{(i,j)}$ as well.

Lemma 6.5 follows from properties about functions having affine-exponential functions as their upper or lower bounds, summarised by the next lemma:

Lemma 6.6.

Assume the setting of Lemma 6.5. Then, we have the following.

• •

$\mathbb{E}_{q(\overline{r})}\left[{\exp(l(r))}\right]<\infty$* for every $l\in\mathcal{A}_{i}$ and every $J$-dimensional normal distribution $q$. *

• •

*For some $l,l^{\prime}\in\mathcal{A}_{i}$, $|\mathit{poly}(\overline{r})|\leq\exp(l(r))$ and $|\mathit{{nn}}(\overline{r})|\leq l^{\prime}(r)$ for all $r\in A_{i}$. *

• •

For some affine functions $l,l^{\prime}:\mathbb{R}\to\mathbb{R}$, $\exp(l(|v|))\leq\mathrm{softplus}(v)\leq\exp(l^{\prime}(|v|))$ for all $v\in\mathbb{R}$.

• •

For every $l_{1},l_{2}\in\mathcal{A}_{i}$ and every $\odot\in\{+,-,\times,/,\max\}$, there exists $l\in\mathcal{A}_{i}$ such that $\exp(l_{1}(r))\odot\exp(l_{2}(r))\leq\exp(l(r))$ for all $r\in A_{i}$.

6.3. Condition for the requirements R4-R6

Assume that the model $C$ and the guide $D_{\theta}$ satisfy our assumption and condition in the previous two subsections. Our condition for the requirements R4-R6 is given below:

[TABLE]

Theorem 6.7.

If both our assumption in §6.1 and the condition (13) hold, then the condition (17) implies the requirements R4-R6.

The proof of the theorem uses the following nontrivial result (Klenke, 2014, Theorem 6.28) about exchanging differentiation and integration, a consequence of the dominated convergence theorem.

Theorem 6.8.

Let $V\subset\mathbb{R}$ be an open interval, and $(X,\Sigma,\mu)$ be a measure space. Suppose that a measurable function $f:V\times X\to\mathbb{R}$ satisfies the following conditions: (a) for all $v\in V$, the integral $\int\mu(\mathrm{d}x)\,f_{v}(x)$ is well-defined; (b) for almost all $x\in X$ (w.r.t. $\mu$) and all $v\in V$, the partial derivative $\nabla_{v}f_{v}(x)$ with respect to $v$ is well-defined;121212A more popular notation is $(\partial f_{v}(x))/(\partial v)$, but we opt for $\nabla_{v}f_{v}(x)$ to avoid clutter. (c) there is a measurable function $h:X\to\mathbb{R}$ such that $\int\mu(\mathrm{d}x)\,h(x)$ is well-defined and $\left|\nabla_{v}f_{v}(x)\right|\leq h(x)$ for all $v\in V$ and almost all $x\in X$ (w.r.t. $\mu$). Then, for all $v\in V$, both sides of the below equation are well-defined, and the equality holds:

[TABLE]

Note that the theorem ensures not only the validity of interchanging differentiation and integration, but also the differentiability of $\int\mu(\mathrm{d}x)\,f_{v}(x)$ (w.r.t. $v$) and the integrability of $\nabla_{v}f_{v}(x)$ over $x\in X$.

Our condition in (17) is sufficient but not necessary for the requirements R4-R6 to hold, in particular for the objective in (5) to have well-defined partial derivatives in $\theta$. However, its violation is a good indication of a potential problem. The following example illustrates this point.

Example 6.9.

Consider a model $C$ and a guide $D_{\theta}$ defined as follows:

[TABLE]

where $E[\theta]\equiv\mathsf{{relu}}(\theta)+2$ and ${\rm relu}(v)\triangleq{\mathds{1}_{[{v\geq 0}]}}\cdot v$. Such $E[\theta]$ can definitely appear in machine learning applications, once a guide starts to use neural networks with parameters $\theta$. Let $A=[\{``a"\}\to\mathbb{R}]$ and $\sigma(\theta)\triangleq{\rm relu}(\theta)+2$. Then, $\mathit{{dens}}(C,s_{I})(r)={\mathds{1}_{[{r\in A}]}}\cdot\mathcal{N}(r(``a");0,1)$ and $\mathit{{dens}}(D_{\theta},s_{I})(r)={\mathds{1}_{[{r\in A}]}}\cdot\mathcal{N}(r(``a");0,\sigma(\theta)).$ Note (17) is violated: $\sigma$ is non-differentiable at $\theta=0$. A simple calculation shows:

[TABLE]

Hence, the objective in (5) does not have a well-defined partial derivative in $\theta$ at $\theta=0$. ∎

7. Analysis

In this section, we describe a recipe for building a static analysis that automatically discharges some of the assumptions and conditions given in §6. The recipe ensures that the constructed static analyses are sound with respect to the density semantics in §4.3. We illustrate it by describing four static analyses for verifying the model-guide support match (the requirement R1), the guide-parameter differentiability (the requirement R2), the condition (13), and the condition (17). The analysis for the model-guide support match has been developed significantly more for the Pyro programming language, and applied to analyse realistic examples of the language. This fully-blown analysis and our experiments will be described in §8.

Throughout this section, we assume that the parameters $\theta$ of a guide $D_{\theta}$ are included in $\mathit{{Var}}$, and are only read by $D_{\theta}$ and not accessed by a model $C$. What we used to call $s_{I}$ will be the part of the store for variables in $\mathit{{Var}}\setminus\theta$, and what we used to write $\theta$ will correspond to the other part for $\theta$.

7.1. A generic program analysis framework

Our recipe is for building a static analysis that infers information about the state transformation of a given command. It is similar to the conventional methodology for building a so called relational static analysis, which also attempts to find information about the relationship between input and output states of a given command. However, our recipe diverges from the convention in one important point: while the abstract states of conventional relational analyses represent relations on states, we let abstract states directly express sets of concrete state transformers. This departure from the convention is due to the difficulty of using relations for expressing properties of state transformers that we desire. For instance, we could not express a set of functions with a certain type of differentiability using relations.

Recall the domain $\mathcal{D}$ in (3), and the notion of admissible subset from domain theory: $D_{0}\subseteq\mathcal{D}$ is admissible if it contains $\bot$ and is closed under taking the limits of $\omega$-chains in $D_{0}$.

Our recipe assumes an abstraction instance defined by the following items:

• •

An abstract domain, i.e., a set $\mathcal{T}^{\sharp}$ with a designated element $\bot^{\sharp}$.

• •

A concretisation function, i.e., a function $\gamma:\mathcal{T}^{\sharp}\to\mathcal{P}(\mathcal{D})$ such that for every $t\in\mathcal{T}^{\sharp}$, $\gamma(t)$ is an admissible subset of $\mathcal{D}$. Note that the concretisation interprets each abstract element $t$ as a set of concrete transformers in $\mathcal{D}$. The admissibility is imposed to enable the sound analysis of loops.

• •

A widening operator $\mathit{{widen}}:\mathcal{T}^{\sharp}\times\mathcal{T}^{\sharp}\to\mathcal{T}^{\sharp},$ such that for all $t_{1},t_{2}\in\mathcal{T}^{\sharp}$ and $i\in[2]$, $\gamma(t_{i})\subseteq\gamma(\mathit{{widen}}(t_{1},t_{2}))$ and for every sequence $\{t_{n}\}_{n\geq 1}$ in $\mathcal{T}^{\sharp}$, its widened sequence $\{t^{\prime}_{n}\}_{n\geq 1}$, defined by $t^{\prime}_{1}\triangleq t_{1}$ and $t^{\prime}_{n+1}\triangleq\mathit{{widen}}(t^{\prime}_{n},t_{n+1})$ for $n\geq 1$, has an index $m$ such that $t^{\prime}_{m}=t^{\prime}_{m+1}$.

• •

An abstract conditional operator for every expression $E$, that is, a function $\mathit{{cond}}(E)^{\sharp}:\mathcal{T}^{\sharp}\times\mathcal{T}^{\sharp}\to\mathcal{T}^{\sharp}$ such that for all $t_{1},t_{2}\in\mathcal{T}^{\sharp}$ and $g_{1},g_{2}\in\mathcal{D}$, if $g_{1}\in\gamma(t_{1})$ and $g_{2}\in\gamma(t_{2})$, then \big{(}\lambda(s,r).(\text{if}\ ({\llbracket E\rrbracket}s{=}\mathit{{true}})\ \text{then}\ g_{1}(s,r)\ \text{else}\ g_{2}(s,r))\big{)}\in\gamma(\mathit{{cond}}(E)^{\sharp}(t_{1},t_{2})).

• •

An abstract composition operator $\circ^{\sharp}:\mathcal{T}^{\sharp}\times\mathcal{T}^{\sharp}\to\mathcal{T}^{\sharp}$ such that for all $t_{1},t_{2}\in\mathcal{T}^{\sharp}$ and $g_{1},g_{2}\in\mathcal{D}$, if $g_{1}\in\gamma(t_{1})$ and $g_{2}\in\gamma(t_{2})$, then $g_{2}^{\ddagger}\circ g_{1}\in\gamma(t_{2}\circ^{\sharp}t_{1}).$

• •

For all expressions $E_{0},E_{1},E_{2}$ and for all variables $x$, the abstract elements $\mathit{{skip}}^{\sharp}$, $\mathit{{update}}(x,E_{0})^{\sharp}$, $\mathit{{sample}}(x,S,E_{1},E_{2})^{\sharp}$, and $\mathit{{score}}(E_{0},E_{1},E_{2})^{\sharp}\,\in\,\mathcal{T}^{\sharp}$ such that

[TABLE]

Given these data, we define the static analysis ${\llbracket C\rrbracket}^{\sharp}\in\mathcal{T}^{\sharp}$ of a command $C$ in Figure 9.

Here the $(\mathit{{wfix}}\ T)$ is the usual widened fixed point $t_{\mathit{fix}}$ of $T$, which is defined as the first element $t_{m}$ with $t_{m}=t_{m+1}$ in the widened sequence $(t_{n})_{n\geq 1}$ where $t_{1}\triangleq\bot^{\sharp}$ and $t_{n+1}\triangleq\mathit{{widen}}(t_{n},T(t_{n}))$.

Theorem 7.1 (Soundness).

For all commands $C$, we have ${\llbracket C\rrbracket}_{d}\in\gamma({\llbracket C\rrbracket}^{\sharp})$.

In the rest of this section, we instantiate this framework into four static analysis instances. In each case, we describe the abstract domain, the abstract bottom element, and the concretisation function. Moreover, in the first two cases, we detail the transfer functions. In the following, for a tuple $(s^{\prime},r^{\prime},w^{\prime},p^{\prime})\in\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty)$, we use the subscripts $-_{s}$, $-_{r}$, $-_{w}$, and $-_{p}$ to denote its components. For instance, $(s^{\prime},r^{\prime},w^{\prime},p^{\prime})_{s}=s^{\prime}$ and $(s^{\prime},r^{\prime},w^{\prime},p^{\prime})_{r}=r^{\prime}$.

7.2. Analysis for the model-guide match

The first instance analysis finds information about the names of sampled random variables. Such information can be used for discharging the requirement R1, the correspondence between the support of a model and that of a guide. The analysis is based on the below abstraction:

[TABLE]

where $[]$ denotes the empty random database. A typical abstract element in $\mathcal{T}^{\sharp}$ is a set of names $K$, which represents concrete commands sampling random variables in $K$. The domain $\mathcal{T}^{\sharp}$ contains $\bot^{\sharp}$ and $\top^{\sharp}$ to express two extreme cases, the set containing only one command that always returns $\bot$, and the set of all commands.

Sound abstract operations can be derived from the density semantics and from the abstraction following the abstract interpretation methodology (Cousot and Cousot, 1977):

[TABLE]

We can use the resulting analysis to discharge the requirement R1. We just need to run it on both $C$ and $D_{\theta}$, and check whether ${\llbracket C\rrbracket}^{\sharp}={\llbracket D_{\theta}\rrbracket}^{\sharp}=K$ for some $K\in\mathcal{P}(\mathit{{Str}})$. The positive answer implies the requirement R1, because all the random variables are drawn from the normal distribution. Our extension of this analysis for Pyro (§8) does not rely on this exclusive use of the normal distribution, and tracks information about the type of distribution of each random variable and state properties, so as to prove the model-guide support match for realistic Pyro programs.

7.3. Analysis for the guide parameter differentiability

The second instance analysis aims at proving the differentiability of the density of a guide $D_{\theta}$ with respect to its parameters $\theta$. It infers the continuous partial differentiability of multiple functions with respect to variables in the input state. The analysis is defined by the below abstraction:

[TABLE]

By “$C^{1}$”, we mean that the relevant function is continuously differentiable. Being a $\mathbb{R}$-valued function in the last part requires that $g(s[z\mapsto v],r)$ be never $\bot$. An $(X,Y,R)$ in $\mathcal{T}^{\sharp}$ expresses a property of a transformer $g\in\mathcal{D}$ (which can be viewed as semantic command) that $g$ does not change variables in $X$, its density is $C^{1}$ with respect to each variable in $Y$ in the input state, and for each $(z,u)\in R$, it assigns a real value to $u$ in the output state in a $C^{1}$ manner with respect to $z$ in the input state.

We now define the abstract operations induced by this abstraction. Given an expression $E$, we let $\mathcal{V}(E)$ denote the set of variables that occur in $E$, and we write $\mathcal{C}^{1}(E)$ for the set of variables with respect to which ${\llbracket E\rrbracket}$ is $C^{1}$ (based on classical differentiability rules). The definitions below follow from general principles such as the multivariate chain rule and account for discontinuities induced by conditions which break differentiability.

[TABLE]

To discharge the differentiability requirement R2, we need to run this analysis on a guide $D_{\theta}$. If the $Y$ component of the analysis result contains all the parameters $\theta$ (i.e., there exists $(X,Y,R)$ such that ${\llbracket D_{\theta}\rrbracket}^{\sharp}=(X,Y,R)$ and $\theta\subseteq Y$), then the requirement R2 is met.

7.4. Analysis for condition (17)

The third analysis extends the second by tracking and checking more properties. Its aim is to prove the condition (17). Just like the second analysis, it infers information about the continuous partial differentiability of multiple functions involving the output state and the density. Also, it checks whether the density of a given command $C_{0}$ has the form

[TABLE]

for some finite $M$, and some $A_{i}$, $K_{i}$, $\mu_{(i,\alpha)}$ and $\sigma_{(i,\alpha)}$, and if so, it tracks properties of the $\mu_{(i,\alpha)}$ and $\sigma_{(i,\alpha)}$. Here is the abstraction for the analysis:

[TABLE]

The abstract operations are similar to those for the differentiability analysis. Thus, we omit their definitions. We can use the analysis to prove the condition (17). We just need to run the analysis on a guide $D_{\theta}$ and check whether ${\llbracket D_{\theta}\rrbracket}^{\sharp}$ is not $\top^{\sharp}$. If so, the condition holds.

7.5. Analysis for condition (13)

The last instance is a static analysis that aims at proving the condition (13). The analysis checks whether a given command $C_{0}$ has a density of the following form:

[TABLE]

for some finite $M$, and some $A_{i}$, $K_{i}$, $N_{i}$, $\mu^{\prime}_{(i,\alpha)}$, $\sigma^{\prime}_{(i,\alpha)}$, $\mu^{\prime\prime}_{(i,j)}$, and $\sigma^{\prime\prime}_{(i,j)}$. If so, it tracks whether the $\mu^{\prime}_{(i,\alpha)}$, $\sigma^{\prime}_{(i,\alpha)}$, $\mu^{\prime\prime}_{(i,j)}$, and $\sigma^{\prime\prime}_{(i,j)}$ and some other functions can be bounded by affine exponentials on the input variables. The abstraction for the analysis is as follows:

[TABLE]

Here we use the notation $\{|s(x)|\}_{x\in Y}$ to mean a partial map from variables $x$ in $Y$ to values $|s(x)|$. The abstract operations that derive from this abstraction are quite similar to those of the differentiability analysis, therefore we do not detail them. To verify the condition (13), we run the analysis on a model $C$ and check whether ${\llbracket C\rrbracket}^{\sharp}\neq\top^{\sharp}$. If the check succeeds, the condition holds.

8. A Static Analysis for Model-Guide Pairs in Pyro and its Evaluation

We present a static analysis that can verify the support correspondence for Pyro model-guide pairs. The analysis extends the first instance of the framework presented in §7. Our presentation focuses on the most salient aspects of this analysis and experimental results.

8.1. Some features of Pyro programs

Pyro is a probabilistic programming framework based on Python and PyTorch. It supports a wide spectrum of distributions and neural networks, and features commands for sampling and scoring. It comes with multiple inference engines, including SVI. We chose Pyro over other probabilistic programming languages (e.g., Anglican) because unlike most of other languages, in Pyro, SVI algorithms are considered as main inference engines and neural networks can be used together with probabilistic programming, which leads to more interesting examples.

In Pyro programs, random-variable names are often created at runtime, for instance by concatenating a fixed string with a sequence of dynamically-generated indices.

Example 8.1 (Use of indexed random-variable names).

The code excerpt (of a model or guide program in Pyro) below samples $N\times M$ instances of independent random variables, and names them with $``x\_1\_1",\ldots,``x\_N\_M"$.

⬇

  for i in range(1,N+1):

    for j in range(1,M+1):

      val = pyro.sample("x_{}_{}".format(i,j), Normal(m,d))  ∎

Since Pyro is based on PyTorch and is adapted to implement data-science applications, Pyro programs heavily use multidimensional arrays, called tensors, and operations over them from PyTorch, in particular, element-wise operations and broadcasting. More precisely, when a binary operation (such as addition) is applied to two tensors of identical size, it outputs a new tensor of the same size where each output element is computed separately from the corresponding input elements. This principle often makes it possible to parallelise computations. Broadcasting occurs when a binary operation is applied to two tensors of different dimensions that can somehow be unified. Intuitively, it introduces and duplicates dimensions to produce tensors of identical dimensions, so that binary operations can be performed element-wise.

Tensors are also heavily used for sampling, which has several consequences. First, it means that an analysis targeted at Pyro programs should be able to track information about at least the dimensions of sampled tensors. Second, the dimensions of sampled tensors are grouped such that each of these groups has a different property with respect to probabilistic independence of tensor components. Pyro inference engines exploit this property for optimisation, for instance, via subsampling, but for this optimisation to be valid, the grouping of dimensions in a model should match that of a guide. Our analysis tracks information about dimension grouping of sampled tensors and the use of the Pyro construct called plate, which enables the optimisation just mentioned.

8.2. Abstract domain

We extend our analysis from §7 so that it tracks not just how each part of a given program transforms states, but also which states can reach at each program point. We need the latter to get information about random variables with dynamically generated names that is precise enough for our verification task. To describe states, we rely on an abstract domain that consists of multiple subdomains (combined by product). Among them, the key part is $\mathit{{RDB}}^{\sharp}\triangleq[\mathit{{Str}}\rightharpoonup_{\mathit{{fin}}}\{\top\}\cup(\mathcal{P}_{\mathit{{fin}}}(\mathit{{Zone}}^{\sharp})\times\mathit{{Dist}}^{\sharp})]$, where $\rightharpoonup_{\mathit{{fin}}}$ denotes a finite partial map. $\mathit{{Dist}}^{\sharp}$ is an abstract domain whose element represents a set of elementary probability distributions, such as the standard normal distribution. The other $\mathit{{Zone}}^{\sharp}$ is our custom domain for zones, which express higher-dimensional rectangles in $\mathbb{N}^{n}$.

An element $r^{\sharp}\in\mathit{{RDB}}^{\sharp}$ means a set of concrete random databases $r$ for each concrete store $s$. The domain of $r$ consists of names that are obtained by concatenating baseline strings $\alpha$ in $\mathrm{dom}(r^{\sharp})$ with index sequences. If $r^{\sharp}$ maps a baseline string $\alpha$ to $\top$, it does not give any information about index sequence and a probability distribution used for the random variable. Otherwise, the second component of $r^{\sharp}$ specifies the information about the distribution, and the first component a region of $\mathbb{N}^{n}$ that contains the index sequences used in the names. The region is described by a finite subset $\{Z_{1}^{\sharp},\ldots,Z_{n}^{\sharp}\}$ of $\mathit{{Zone}}^{\sharp}$, which means the disjoint union of the rectangles represented by the $Z_{i}^{\sharp}$. To emphasise this union interpretation, we write $Z_{1}^{\sharp}\cup\ldots\cup Z_{n}^{\sharp}$ for the subset. The following table summarises the definitions of our abstract domain, based on $\mathit{{Zone}}^{\sharp}$:

[TABLE]

A higher-dimensional rectangular zone $Z^{\sharp}$ is described by a finite sequence of intervals $I^{\sharp}$, each of which is made of two bounds. A bound $B^{\sharp}$ may be defined as one or two constraints which express that this bound is equal to a constant, or to a variable plus a constant, or both. This intuitive denotation defines a concretisation function $\gamma_{r}:\mathit{{RDB}}^{\sharp}\to\mathcal{P}(\mathit{{State}})$.

Example 8.2.

Consider the code of Example 8.1. After $i-1$ iterations in the main loop and $j$ iterations in the last execution of the inner loop, we expect the analysis to come up with the invariant, $[``x"\mapsto([1,i{-}1]\times[1,M]\cup[i,i]\times[1,j],\,\mathrm{normal}(m,d))]$. In turn, at the exit of the main loop, we expect the analysis to infer the invariant, $[``x"\mapsto([1,i{=}N]\times[1,j{=}M],\,\mathrm{normal}(m,d))]$. ∎

In addition to these constraints over random databases, our analyser also uses an abstraction that maintains typing information and numerical constraints over variables. We do not fully formalise these constraints as the overall structure of the domain relies on a classical reduced product.

Finally, we describe the combination of the above state abstraction with the abstraction of §7.2. More precisely, we start with the abstract domain exposed in §7.2 (which we denote by $\mathcal{T}^{\sharp}$) and build a novel abstract domain $\mathcal{T}^{\sharp}_{s}$ that also satisfies the requirements of §7.1. We let the set of abstract elements be $\mathcal{T}^{\sharp}_{s}\triangleq[RDB^{\sharp}\rightarrow\mathcal{T}^{\sharp}\times RDB^{\sharp}]$. Intuitively, such an element maps an input abstract random database into a set of functions together with an over-approximation of their output, when applied to this abstract input. The concretisation below formalises this: for all $t_{s}\in\mathcal{T}^{\sharp}_{s}$,

[TABLE]

8.3. Computation of loop invariants

Although our static analysis for Pyro requires a state abstraction, its principles and structure are similar to those of the general analysis shown in §7. In the following, we first describe the integration of the state abstraction in the analysis of §7.2.

The abstract operations in $\mathcal{T}^{\sharp}_{s}$ can all be derived by lifting those in $\mathcal{T}^{\sharp}$ into functions. We illustrate this for the abstract composition $\circ^{\sharp}_{s}$ for $\mathcal{T}^{\sharp}_{s}$. Recall the operator $\circ^{\sharp}$ for $\mathcal{T}^{\sharp}$ in §7.2. Given $t_{s},t^{\prime}_{s}\in\mathcal{T}^{\sharp}_{s}$,

[TABLE]

The other abstract operators over $\mathcal{T}^{\sharp}_{s}$ are defined in a similar manner.

As in most static analysis problems, the computation of precise loop invariants requires a carefully designed widening operator. In the case of zones, the analysis needs to generalise information about the bounds. We assume $r^{\sharp}_{0},r^{\sharp}_{1}\in\mathit{{RDB}}^{\sharp}$ are abstract random databases, and present the main steps in the computation of the widening $\mathit{{widen}}(r^{\sharp}_{0},r^{\sharp}_{1})$.

1. (i)

For each $r^{\sharp}_{i}$, we fuse zones $Z^{\sharp}_{0},Z^{\sharp}_{1}$ in $r^{\sharp}_{i}$ into a single zone when $Z^{\sharp}_{0}$ and $Z^{\sharp}_{1}$ are the same except for one component and they represent adjacent high-dimensional rectangles. The following rewriting illustrates this step.

[TABLE] 2. (ii)

We generalise the intervals of corresponding zones in $r^{\sharp}_{0}$ and $r^{\sharp}_{1}$ using a weak version of unification on their bounds. The bounds of corresponding intervals in two such zones survive this step in a weakened form if they share at least one syntactically equal expression. Otherwise, the bounds are dropped, which causes the introduction of $\top$ into abstract random databases. The following rewriting instances illustrate this bound generalisation.

[TABLE]

The analysis also applies standard widening techniques to typing information and numerical constraints mentioned in §8.2. Finally, a model-guide pair can be verified if and only if their analyses return equivalent abstract random databases, without any name mapped into $\top$.

8.4. Experimental evaluation

We have implemented a prototype analyser and carried out experiments so as to assess the effectiveness of our analysis to verify the support correspondence for Pyro model-guide pairs.131313Code is available at https://github.com/wonyeol/static-analysis-for-support-match. More precisely, we evaluated the following three research questions:

1. (i)

Can our analysis discover incorrect model-guide pairs in realistic probabilistic programs? 2. (ii)

Can the analysis verify correct model-guide pairs despite the complexity of the definition? 3. (iii)

Is the analysis efficient enough so that it can cope with realistic probabilistic programs?

Benchmark programs.

We took example programs from the Pyro webpage that can be handled by standard SVI engines or can be naturally rewritten to be so. Those engines do not use Pyro’s recent vectorised enumeration-based optimisation (enabled by the option TraceEnum_ELBO), and are run with the option Trace_ELBO. We made this choice because the optimisation is newly introduced and whether it is used or not changes the set of random variables used in a Pyro program.

We have applied our analysis to two groups of Pyro programs. The first is the Pyro regression test suite (Uber AI Labs, 2019b), which comprises 66 examples exercising the basic programming patterns expected in Pyro programs. While each of these is small, they cover many standard patterns of both correct and incorrect model-guide pairs. Among these, 33 test cases use only TraceEnum_ELBO and fall outside the scope of our analysis, and 6 of the remaining 33 test cases come with two subsampling options. Hence, we can consider 39 experiments based on this regression suite.

The second is a collection of examples from the Pyro webpage (Uber AI Labs, 2019a). The webpage features 18 Pyro examples, among which 9 involve model-guide pairs (other examples use automatic guide generation, or do not perform SVI at all). Out of these 9 examples, 5 use Trace_ELBO, and three can be converted naturally into Trace_ELBO versions although they do not use Trace_ELBO. Thus, 8 Pyro examples fall within the scope of our analysis. These examples correspond to advanced probabilistic models141414The models include variational autoencoder (VAE) (Kingma and Welling, 2014), semi-supervised VAE (Kingma et al., 2014), deep exponential family (Ranganath et al., 2015), attend-infer-repeat (Eslami et al., 2016), deep Markov model (Krishnan et al., 2017), inference compilation (Le et al., 2017), and amortised latent Dirichlet allocation (Srivastava and Sutton, 2017).

from the machine-learning literature, most of which use sophisticated neural networks and probabilistic modelling. Table 1 describes the structure of these 8 examples in detail.

Prototype analyser and results.

Our analyser is implemented in OCaml, and supports the main data-structures and operations defined by Python, PyTorch, and Pyro. In particular, it precisely abstracts the shape of PyTorch tensor objects, the shape transformation information of PyTorch neural-network-related objects, the automatic broadcasting information of Pyro plate objects, and the shape of allocated indices for sample names, using the zone abstraction described above. It also supports common Pyro probability distributions, and can precisely cope with standard Pyro and PyTorch operations manipulating the Pyro distribution objects and PyTorch tensor objects. While our prototype supports a wide range of Python, PyTorch, and Pyro features, we point out that we did not implement a static analysis for the full Python (plus PyTorch and Pyro) language (e.g., no support for classes, dictionaries, named tuples, and user-defined functions).

The analysis results are summarised in Table 2 and are discussed in detail in the following. Run-times were measured on an Intel Core i7-7700 machine running Linux Ubuntu 16.04.

Discovery of invalid model-guide pairs.

The analysis rejected two Pyro examples, br and lda, as incorrect due to an invalid model-guide pair. br is the Bayesian regression example discussed in §2.

For br, the analysis discovers that a random variable sigma is sampled from Uniform(0.,10.) in the model $C$, but from Normal(..) in the guide $D_{\theta}$ (Figure 2(a)). Since the support of sigma in $D_{\theta}$ is not a subset of that in $C$ (i.e., $\mathbb{R}\not\subseteq[0,10]$), the requirement R1 is violated. Thus, the SVI objective, $\mathrm{KL}(D_{\theta}\|C)$, is undefined, and br has an invalid model-guide pair.

For lda, the analysis discovers that a random variable doc_topics is sampled from Dirichlet(..) in the model $C$, but from Delta(..) in the guide $D_{\theta}$. Since the reference measures of doc_topics in $C$ and $D_{\theta}$ are different (the Lebesgue measure vs. the counting measure), $\mathrm{KL}(D_{\theta}\|C)$ cannot be computed by (4). For this reason, lda has an invalid model-guide pair. Our analyser tracks the reference measure implicitly by regarding the support of any distribution with Lebesgue measure, as disjoint from that of any distribution with counting measure (which is sound due to the aforementioned reason), and this allowed us to discover the issue of lda.

In both cases, the found correctness issues are unknown before and subtle. In particular, it turned out that lda, though incorrect when viewed as an example of SVI, is a valid implementation, because it performs not SVI but variational expectation-maximisation (EM) (Neal and Hinton, 1998).151515We thank an anonymous reviewer and Eli Bingham for pointing this out and explaining it in detail. The lda uses an SVI engine just to solve an optimisation problem in variational EM (not to do SVI), and uses the Delta distribution to perform a particular realisation of the M-step in variational EM.

Verification of probabilistic programs relying on model-guide pairs.

Among the Pyro test suite, the analysis successfully verifies 31 examples among 39. Interestingly, two of these 31 successful validations, highlighted in Table 2(a), correspond to cases that were flagged as “invalid model-guide pairs” in the Pyro git repository. Upon inspection, these two examples turn out to be correct.

On the other hand, 8 examples from the Pyro test suite could not be verified due to the crashes of the analyser. One of these failures is due to the need to reason more precisely about the content of a for loop (e.g., using some partitioning techniques), and seven are due to the use of plates with subsampling, as ranges for for loops. Therefore, these failures could be resolved using existing static analysis techniques and a more precise handling of the semantics of Python constructions.

Moreover, 6 Pyro examples (among the 8 that we considered) were verified successfully, which means all correct Pyro examples were verified. Finally, we corrected the two examples that were rejected due to invalid model-guide pairs, and these two examples were also successfully verified.

Analysis efficiency.

The analysis returned within a second on each program in the Pyro test suite, and on most of the Pyro examples. The slowest analysis was observed on air, which was analysed within 5 seconds. Most of the Pyro examples sample from and score with distributions of very high dimension arranged in complex tensors, using nested for and plate’s. While they are not large, they present a high degree of logical complexity, that is representative of realistic probabilistic programs. The fact that such programs get analysed within seconds shows that the analysis and the underlying abstract domain to describe zones, sampled dimensions, and distributions can generalise predicates quickly so that precise loop invariants can be computed.

9. Related Work, Conclusion and Limitation

Related work.

As far as we know, the idea of using SVI for probabilistic programs first appeared in (Wingate and Weber, 2013). When further insights into how to create generic (sometimes also called black-box) SVI engines were found (Ranganath et al., 2014; Kucukelbir et al., 2015, 2017), the idea was tried for realistic probabilistic programming languages, such as Stan (Kucukelbir et al., 2015) and Anglican (van de Meent et al., 2016), resulting in impressive case studies (Kucukelbir et al., 2015). The major inference engines for deep probabilistic programming languages (Bingham et al., 2019; Tran et al., 2016, 2018; Siddharth et al., 2017) are all based on SVI nowadays. However, we do not know of any prior work that attempts to reveal implicit assumptions made by these SVI engines and to discharge these assumptions manually or automatically, as we did in this paper.

The correctness of a different type of inference engines based on Monte-Carlo methods has been the subject of study in the PL community. Such engines have clearly formulated correctness conditions from Markov chain theory (Metropolis et al., 1953; Hastings, 1970; Green, 1995; Geyer, 2011), such as ergodicity and correct stationarity. Tools from formal semantics have been employed to show that the inference engines satisfy these conditions (Hur et al., 2015; Kiselyov, 2016; Borgström et al., 2016; Scibior et al., 2018). While looking at different algorithms, some of these work consider more expressive languages than what we used in the paper, in particular, those supporting higher-order functions. One interesting direction is to extend our results to such expressive languages using the ideas from these works, especially the operational technique in (Borgström et al., 2016) and the denotational monad-based technique in (Scibior et al., 2018).

The consequence of having random choice in a programming language has been actively investigated by the semantics researchers from the early days (Kozen, 1981; Jones and Plotkin, 1989; Ehrhard et al., 2014; Toronto et al., 2015; Borgström et al., 2016; Staton et al., 2016; Staton, 2017; Smolka et al., 2017; Heunen et al., 2017; Vákár et al., 2019). Our work uses the technique developed in this endeavour, such as Giry monad and denotational formulation of idealised importance sampling (Staton et al., 2016). Also, just as we connected the measure semantics with the density semantics, (Kozen, 1981) and (Borgström et al., 2016) related two semantics with similar flavours, although the considered languages (with or without continuous distribution and score statements) and the style of semantics (operational vs denotational) are different. Our density semantics uses a reference measure built out of Lebesgue measure as in (Bhat et al., 2012, 2013; Hur et al., 2015). This choice is found to cause an issue when the score statement is treated not as a scoring mechanism but in terms of conditional expectation (Wu et al., 2018). How to resolve this matter is still open.

Static analyses for probabilistic programming languages or languages with probabilistic choice typically attempt to check probabilistic properties (Monniaux, 2000; Cousot and Monerau, 2012; Chakarov and Sankaranarayanan, 2013; Wang et al., 2018), or perform posterior inference (Chaganty et al., 2013; Wang et al., 2018), or find information useful for posterior inference (Nori et al., 2014). For instance, (Monniaux, 2000) defines a probabilistic abstract interpretation framework, which is applied to estimate the failure probability (Monniaux, 2001). More recently, (Cousot and Monerau, 2012) sets up a general framework to design probabilistic abstract interpretations, which lift conventional static analysis tools to a probabilistic setup, and (Wang et al., 2018) proposes a framework for analysing probabilistic programs based on hyper-graphs with probabilistic actions on edges. Our program analyses aim at a different type of verification tasks, namely, proving the safety requirements imposed by the SVI engines. Static analyses for checking the continuity properties of programs are proposed in (Chaudhuri et al., 2010). Some techniques used in those analyses may help track the kind of smoothness properties considered in this paper precisely.

Conclusion.

In this paper, we have set up a semantic model of probabilistic programming languages that allows to reason about the correctness of subtle optimisation algorithms such as SVI. In particular, we have identified a set of conditions that guarantee that these algorithms will produce meaningful results. These conditions simplify the verification of correctness, and open the way towards the automation of this verification. To evaluate the effectiveness of our framework, we have completed the design of a static analysis to validate the correspondence assumptions for the model-guide pairs. Our analysis could verify nontrivial implementations as well as uncover nontrivial bugs. Such bugs would lead to incorrect results, which would be very difficult to diagnose.

Limitation.

We described four static analyses (§7.2-§7.5) for discharging the requirements R1-R6, and developed one of the four analyses further (§8.2-§8.4) to build a prototype analyser for actual Pyro programs. Here we clarify the underlying assumptions and limitations of these analyses.

The static analysis for model-guide support match (§7.2) was implemented into our static analyser for Pyro (§8.4) but with additional extensions (§8.2-§8.3), so that it does not make the assumption in §6.1; the assumption was introduced mainly to develop analyses for the requirements R3-R6. Hence, our static analyser handles both continuous and discrete random variables.

On the other hand, other analyses (§7.3-§7.5) are unimplemented and make the assumption in §6.1. We expect that the assumption can be relaxed, without much difficulty, to permit continuous or discrete distributions having finite entropy and moments of all degrees, because our proofs of theorems and lemmas bound various expectations mostly by entropies and moments. It would be more challenging, however, to relax the assumption to allow distributions not having finite entropy and moments of all degrees, or models having unboundedly many random variables (of any kind). In particular, addressing the later limitation might require techniques developed for proving that probabilistic systems has finite expected execution time.

This paper focuses on a particular optimisation objective (5) for SVI, but we point out that several other optimisation objectives have been proposed for different inference or learning algorithms, such as variational EM (e.g., lda in §8.4) and importance weighted autoencoder (IWAE) (Burda et al., 2016). One interesting research direction is to develop techniques for verifying the well-definedness of these optimisation objectives.

Acknowledgements.

We thank Fritz Obermeyer and Eli Bingham for explaining the subtleties of Pyro and suggesting us to try the Pyro regression test suite. Sam Staton, Ohad Kammar and Matthijs Vákár helped us to understand the semantics of recursion in the probabilistic programming languages better. Lee, Yang and Yu were supported by the Engineering Research Center Program through the National Research Foundation of Korea (NRF) funded by the Korean Government MSIT (NRF-2018R1A5A1059921), and also by Next-Generation Information Computing Development Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT (2017M3C4A7068177).

Appendix A Proofs of Lemmas and Theorems in §4

Lemma 4.1. For all expressions $E$, $B$, and $S$, their semantics ${\llbracket E\rrbracket}$, ${\llbracket B\rrbracket}$ and ${\llbracket S\rrbracket}$ are measurable functions from $\mathit{{Store}}$ to $\mathbb{R}$, $\mathbb{B}$ and $\mathit{{Str}}$, respectively.

Proof.

We prove the lemma by structural induction. We present the proof for $E$ only, because those for $B$ and $S$ are similar. When $E$ is a constant, ${\llbracket E\rrbracket}$ is a constant function and so it is measurable. When $E$ is a variable, ${\llbracket x\rrbracket}$ is a projection of the $x$ component of the given store $s$. Since $\mathit{{Store}}$ is given the product $\sigma$-algebra, the projection ${\llbracket x\rrbracket}$ is measurable. Finally, assume that $E$ is $f(E_{0},\ldots,E_{n-1})$. Then, by induction hypothesis, the ${\llbracket E_{i}\rrbracket}$ are measurable functions from $\mathit{{Store}}$ to $\mathbb{R}$. Thus,

[TABLE]

is measurable. Furthermore, ${\llbracket f\rrbracket}:\mathbb{R}^{n}\to\mathbb{R}$ is measurable. Thus, their composition is measurable as well. This implies the measurability of ${\llbracket E\rrbracket}$ because

[TABLE]

∎

Lemma 4.2. $(\mathcal{K},\sqsubseteq)$* is an $\omega$-complete partial order with the least element $\bot$.*

Proof.

The fact that $\sqsubseteq$ is a partial order immediately follows from its definition. The least element is the function $\lambda\sigma.\,\lambda A.\,0$. It remains to show the $\omega$-completeness property. Consider an $\omega$-chain

[TABLE]

Let $\kappa=\sup_{n}\kappa_{n}$. We will show that $\kappa$ is the desired upper bound.

Pick a state $\sigma$. Then, $\kappa(\sigma)(A)$ is non-negative for all measurable subsets $A$, and $\kappa(\sigma)(\emptyset)=0$. Furthermore, for every countable collection $\{A_{j}\}_{j}$ of disjoint measurable subsets of $\mathit{{State}}\times[0,\infty)$,

[TABLE]

The exchange between $\sup$ and $\lim$ in the derivation uses the fact that the sequence involved is increasing. The third equality follows from the countable additivity of $\kappa_{n}(\sigma)$, and the fourth from the monotone convergence theorem applied to the counting measure on the set of natural numbers. Finally,

[TABLE]

We have just shown that $\kappa(\sigma)$ is a subprobability measure.

To show the measurability of $\kappa$, we use the fact that a map $f$ from a measurable space $(X,\Sigma)$ to $\mathrm{Sp}(\mathit{{State}}\times[0,\infty))$ is measurable if and only if for all measurable subsets $A$ of $\mathit{{State}}\times[0,\infty)$ and reals $r\in\mathbb{R}$,

[TABLE]

Pick a measurable subset $A$ of $\mathit{{State}}\times[0,\infty)$ and a real $r\in\mathbb{R}$. Then,

[TABLE]

But the intersection at the end of this derivation gives a measurable set because each $\{\sigma\,\mid\,\kappa_{n}(\sigma)(A)\}$ is measurable and the countable intersection preserves measurability. ∎

Lemma A.1.

Let $\{\kappa_{n}\}_{n}$ be an $\omega$-chain in $\mathcal{K}$. For all non-negative measurable functions $h:\mathit{{State}}\times[0,\infty)\to[0,\infty)$ and all states $\sigma$,

[TABLE]

Proof.

We prove it by a simpler version of the so called monotone class theorem. This version says that in order to show the claimed equality for all non-negative measurable $h$, it suffices to show the following properties:

• •

The claimed equality holds if $h$ is the indicator function $\lambda(\sigma,w).\,{\mathds{1}_{[{(\sigma,w)\in A}]}}$ for a measurable subset $A$ of $\mathit{{State}}\times[0,\infty)$.

• •

If non-negative measurable $h_{1}$ and $h_{2}$ satisfy the claimed equality and $r_{1},r_{2}$ are non-negative real numbers, the linear combination $r_{1}\cdot h_{1}+r_{2}\cdot h_{2}$ also satisfies the equality.

• •

If $\{h_{m}\}_{m}$ is a countable family of non-negative measurable functions such that each $h_{m}$ satisfies the claimed equality and $h_{m}\leq h_{m+1}$ for every $m$, then $\sup_{m}h_{m}$ also satisfies the claimed equality.

We will show these properties one by one.

Let us start with the proof of the first property. Pick a measurable subset $A$ of $\mathit{{State}}\times[0,\infty)$. Then,

[TABLE]

The first and the last equalities use the definition of integration, the second the characterisation of the least upper bound in $\mathcal{K}$, and the third the fact that $\{\kappa_{n}(\sigma)(A)\}_{n}$ is increasing.

For the second property, assume that the claimed equality holds for non-negative measurable functions $h_{1}$ and $h_{2}$. Let $r_{1}$ and $r_{2}$ be non-negative real numbers. Then,

[TABLE]

The first and the last equalities use the linearity of integration, the second the assumption on $h_{1}$ and $h_{2}$, and the third the continuity of constant multiplication and addition.

For the final third property, consider an increasing sequence $\{h_{m}\}_{m}$ in the property. Then,

[TABLE]

The third equality uses the monotone convergence theorem, and the fourth the assumption on the $h_{m}$. ∎

Theorem 4.3. For every command $C$, its interpretation ${\llbracket C\rrbracket}$ is well-defined and it belongs to $\mathcal{K}$.

Proof.

We prove the theorem by induction on the structure of $C$. In the proof, we use concepts and results from the measure theory that go beyond what we reviewed in §3.

We start with the case $C\equiv(x:=E)$. For every $\sigma\in\mathit{{State}}$, ${\llbracket x:=E\rrbracket}(\sigma)$ is the Dirac measure at $(\sigma[x\mapsto{\llbracket E\rrbracket}\sigma_{s}],1)$. Thus, it suffices to show the measurability of the following function $F$ for every measurable subset $A$:

[TABLE]

This holds because the functions

[TABLE]

are both measurable.

The case $C\equiv\mathsf{{skip}}$ follows from our argument for the previous case $C\equiv(x:=E)$, because ${\llbracket\mathsf{{skip}}\rrbracket}={\llbracket x:=x\rrbracket}$.

Next we consider sequential composition. Assume that $C\equiv(C_{0};C_{1})$. By induction hypothesis, we know that ${\llbracket C_{0}\rrbracket}$ and ${\llbracket C_{1}\rrbracket}$ are subprobability kernels. Then, for all $\sigma$ and $A$,

[TABLE]

Furthermore, it is not difficult to show that ${\llbracket C_{0};C_{1}\rrbracket}(\sigma)$ is a measure. The measurability of ${\llbracket C_{0};C_{1}\rrbracket}$ follows from three facts. First, it suffices to show that for all measurable subsets $A$, the function

[TABLE]

is measurable. Second, for every measurable $f$ and every subprobability kernel $\kappa$, the function $x\longmapsto\int\kappa(x)(\mathrm{d}y)f(y)$ is measurable. Second, the functions $(\sigma,w)\longmapsto{\llbracket C_{1}\rrbracket}(\sigma)\otimes\delta_{w}$ and $((\sigma^{\prime\prime},w^{\prime\prime}),w)\longmapsto(\sigma^{\prime\prime},w^{\prime}\cdot w^{\prime\prime})$ are measurable. Here $\delta_{w}$ is the Dirac measure at $w$ and $\otimes$ is the operator for product measure.

We move on to the case of $C\equiv(\mathsf{{if}}\ B\ \{C_{0}\}\ \mathsf{{else}}\ \{C_{1}\})$. We first show that the function ${\llbracket C\rrbracket}$ is a well-defined function from $\mathit{{State}}$ to $\mathrm{Sp}(\mathit{{State}}\times[0,\infty))$, that is, ${\llbracket C\rrbracket}(\sigma)$ is a subprobability measure for every $\sigma\in\mathit{{State}}$. The non-negativity and countable additivity of ${\llbracket C\rrbracket}(\sigma)$ follow from the same properties of ${\llbracket C_{0}\rrbracket}(\sigma)$ and ${\llbracket C_{1}\rrbracket}(\sigma)$. Also,

[TABLE]

It remains to show that ${\llbracket C\rrbracket}$ is measurable. This boils down to showing that

[TABLE]

is measurable for all $r\in\mathbb{R}$ and measurable subset $A$ of $\mathit{{State}}\times[0,\infty)$. But this follows easily from the induction hypothesis and the measurability of ${\llbracket B\rrbracket}$.

The next case is $C\equiv(\mathsf{{while}}\,B\,\{C_{0}\})$. Recall Lemma 4.2, which says that the set of subprobability kernels from $\mathit{{State}}$ to $\mathit{{State}}\times[0,\infty)$, denoted by $\mathcal{K}$, is an $\omega$-cpo when given the pointwise order $\sqsubseteq$. Let $F$ be the function on $\mathcal{K}$ used in the semantics of $C$. That is, for all $\kappa\in\mathcal{K}$, $\sigma\in\mathit{{State}}$ and measurable $A\subseteq\mathit{{Store}}\times[0,\infty)$,

[TABLE]

It suffices to show that $F$ is a continuous function on $\mathcal{K}$. The arguments used to handle the if and sequencing cases show that $F(\kappa)$ is a subprobability kernel. The monotonicity of $F$ follows from the fact that the integration is monotone with respect to the measure and its integrand. It remains to show that $F$ preserves the limits of $\omega$-chains. Consider an $\omega$-chain $\{\kappa_{n}\}_{n}$ of subprobability kernels in $\mathcal{K}$. Then, the least upper bound of this $\omega$-chain is $\kappa=\bigsqcup_{n}\kappa_{n}$. The following calculation shows that $F$ preserves this limit.

[TABLE]

The third equality uses the continuity of integration with respect to the integrating measure and the non-negative integrand (Lemma A.1), and the fourth equality uses the monotone convergence theorem.

The next case is $C\equiv(x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2}))$. The requirement that

[TABLE]

is a measure for every $\sigma$ follows immediately from the semantic equation for this sample statement, which says that ${\llbracket x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})\rrbracket}(\sigma)$ is either the constant-[math] measure or the pushforward of a normal distribution with mean ${\llbracket E_{1}\rrbracket}\sigma_{s}$ and standard deviation ${\llbracket E_{2}\rrbracket}\sigma_{s}$ by some measurable function. In fact, the following calculation shows that ${\llbracket x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})\rrbracket}$ is a subprobability measure:

[TABLE]

Also, since all of $\mathcal{N}$, ${\llbracket S\rrbracket}$, ${\llbracket E_{0}\rrbracket}$ and ${\llbracket E_{1}\rrbracket}$ are measurable,

[TABLE]

is a measurable subset of $\mathit{{State}}$ for all $A$ and $r$. This means that ${\llbracket x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})\rrbracket}$ is a measurable function. From what we have proved, it follows that ${\llbracket x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})\rrbracket}$ is a subprobability kernel as required.

The last case is $C\equiv\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})$. Because of the semantic equation for this score statement, ${\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}(\sigma)$ is a measure for all $\sigma$. Furthermore, ${\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}(\sigma)(A)$ is the product of two indicator functions, and so it is at most $1$. That is, ${\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}(\sigma)$ is a subprobability measure. By the semantic equation for the score statement again, the measurability of the normal density and the measurability of $E_{0}$, $E_{1}$ and $E_{2}$, the set $\{\sigma\,\mid\,{\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}(\sigma)(A)\leq r\}$ is measurable for all $A$ and $r$. This measurability implies that ${\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}$ is a measurable function. By putting all these together, we can derive the required fact that ${\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}$ is a subprobability kernel. ∎

Lemma 4.4. For every measurable $h:\mathit{{RDB}}\times\mathit{{RDB}}\times\mathit{{RDB}}\to\mathbb{R}$, the function $(r,r^{\prime})\longmapsto{\mathds{1}_{[{r\#r^{\prime}}]}}\times h(r,r^{\prime},r\uplus r^{\prime})$ from $\mathit{{RDB}}\times\mathit{{RDB}}$ to $\mathbb{R}$ is measurable.

Proof.

Let ${*}:\mathit{{RDB}}\times\mathit{{RDB}}\to\mathit{{RDB}}$ be the following extension of the partial function $\uplus:\mathit{{RDB}}\times\mathit{{RDB}}\to_{\mathit{{partial}}}\mathit{{RDB}}$ to a total function:

[TABLE]

Then, for all $r,r^{\prime}$,

[TABLE]

We will show that $(r,r^{\prime})\longmapsto{\mathds{1}_{[{r\#r^{\prime}}]}}$ and $(r,r^{\prime})\longmapsto r*r^{\prime}$ are measurable. This is enough because the function on the RHS of the above equation (namely, $(r,r^{\prime})\longmapsto{\mathds{1}_{[{r\#r^{\prime}}]}}\cdot h(r,r^{\prime},r*r^{\prime})$) is constructed by applying pairing, composition and multiplication operators on measurable functions and so it is measurable.

Pick a real number $v\in\mathbb{R}$. We should show that $B=\{(r,r^{\prime})\,\mid\,{\mathds{1}_{[{r\#r^{\prime}}]}}\leq v\}$ is measurable. If $v\geq 1$, the set $B$ is $\mathit{{RDB}}\times\mathit{{RDB}}$, and so it is measurable. If $v<0$, then $B$ is the empty set, so that it is measurable as well. If $0\leq v<1$, then

[TABLE]

which is measurable because the RHS of the above equation is the countable union of measurable sets.

Pick a measurable subset $A$ of $\mathit{{RDB}}$. We have to prove that $B=\{(r,r^{\prime})\,\mid\,r*r^{\prime}\in A\}$ is measurable. If $A=\{[]\}$, then $B=\{(r,r^{\prime})\,\mid\,{\mathds{1}_{[{r\#r^{\prime}}]}}\leq 0.5\}\cup\{([],[])\}$, so that it is measurable. If $[]\not\in A$, then

[TABLE]

But $\{r^{\prime\prime}\in A\,\mid\,\mathrm{dom}(r^{\prime\prime})=K\cup K^{\prime}\}$ is measurable, and the set on the RHS of the above equation is the countable union of measurable sets. Thus, $B$ is measurable in this case as well. The remaining case is that $A$ contains $[]$ and some other element. We can prove this case by splitting $A$ into $\{[]\}$ and $A\setminus\{[]\}$ and using what we proved in the previous two cases. ∎

Lemma 4.5. $(\mathcal{D},\sqsubseteq)$* is an $\omega$-complete partial order and has the least element $a\longmapsto\bot$. Thus, every continuous function $G$ on $\mathcal{D}$ has the least fixed point.*

Proof.

All constant functions from a measurable space to another are measurable. Thus, the constant function $a\longmapsto\bot$ is measurable. It is immediate that this constant function is local and smaller than any other functions in $\mathcal{D}$.

Let $g_{0}\sqsubseteq g_{1}\sqsubseteq g_{2}\sqsubseteq\ldots$ be an $\omega$-chain. We show that the following function is the least upper bound of the sequence in $\mathcal{D}$: for all $s,r$,

[TABLE]

By definition, this function is the least upper bound in $\mathcal{D}$ that we are looking for, if it is measurable and local. To check the locality, suppose that $(\bigsqcup_{i}g_{i})(s,r)=(s^{\prime},r^{\prime},w^{\prime},p^{\prime})$. Then, there exists some $i$ such that $g_{i}(s,r)=(s^{\prime},r^{\prime},w^{\prime},p^{\prime})$. By the locality of $g_{i}$, we have that

[TABLE]

Thus, the definition of $\bigsqcup_{i}g_{i}$ implies that the above property holds when we replace $g_{i}$ by $\bigsqcup_{i}g_{i}$. We have just shown that $\bigsqcup_{i}g_{i}$ is local.

It remains to show that $\bigsqcup_{i}g_{i}$ is measurable. Let $A\subseteq\{\bot\}\cup(\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty))$ be a measurable set. We prove that $(\bigsqcup_{i}g_{i})^{-1}(A)$ is a measurable set.

• •

Case $A=\{\bot\}$:

[TABLE]

The intersection of countably many measurable sets are measurable, and $g_{0},g_{1},\ldots$ are measurable functions. Hence we get the measurability of $(\bigsqcup_{i}g_{i})^{-1}(\{\bot\})$.

• •

Case $\bot\not\in A$:

[TABLE]

The first equivalence follows from the definition of the order $\sqsubseteq$. The union of measurable sets are measurable, and $g_{0},g_{1},\ldots$ are measurable functions. Thus, the desired result follows.

• •

Case $A=\{\bot\}\cup A^{\prime}$ with $\bot\not\in A^{\prime}$:

Then $A^{\prime}$ is also measurable as $\{\bot\}$ is and the set subtraction preserves measurability. Moreover $(\bigsqcup_{i}g_{i})^{-1}(A)=(\bigsqcup_{i}g_{i})^{-1}(\{\bot\})\cup(\bigsqcup_{i}g_{i})^{-1}(A^{\prime})$, and the union of countably many measurable sets is measurable, and so we can apply the two cases above.

∎

The $\omega$-cpo $\mathcal{D}$ permits an operator for composing functions in it. We say that a function

[TABLE]

is local if for all $(s,r,w,p),(s^{\prime},r^{\prime},w^{\prime},p^{\prime})\in\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty)$,

[TABLE]

We say that $g$ is strict if $g(\bot)=\bot$. Let

[TABLE]

Define a partial order $\sqsubseteq$ on $\mathcal{D}^{\ddagger}$ by

[TABLE]

Lemma A.2.

$(\mathcal{D}^{\ddagger},\sqsubseteq)$* is an $\omega$-complete partial order and has the least element.*

Proof.

Let

[TABLE]

Essentially by the same argument as the one used in the proof of Lemma 4.5, the set $\mathcal{E}$ with the pointwise order $\sqsubseteq$ is an $\omega$-cpo and has the least element. But $(\mathcal{D}^{\ddagger},\sqsubseteq)$ and $(\mathcal{E},\sqsubseteq)$ are isomorphic as partially order sets. Thus, $(\mathcal{D}^{\ddagger},\sqsubseteq)$ has the property claimed by the lemma. ∎

The following lemma says that the lifting $g^{\ddagger}$ of $g\in\mathcal{D}$ in §4.3 is a well-defined function from $\mathcal{D}$ to this new domain.

Lemma A.3.

The function $g\longmapsto g^{\ddagger}$ from $\mathcal{D}$ to $\mathcal{D}^{\ddagger}$ is well-defined and continuous.

Proof.

We first show that the function is well-defined. Let $g\in\mathcal{D}$. Then, $g^{\ddagger}(\bot)=\bot$ by the definition of $g^{\ddagger}$. The measurability of $g^{\ddagger}$ follows from the measurability of $g$, projection, and multiplication and the fact that both the pairing of measurable functions (into a function to a product measurable space) and the function composition preserve measurability. The locality of $g^{\ddagger}$ follows from that of $g$.

Next, we show the claimed continuity. The monotonicity is an immediate consequence of the definitions of the $-^{\ddagger}$ operator and the orders on $\mathcal{D}$ and $\mathcal{D}^{\ddagger}$. It remains to show the limit preservation. Let $g_{0}\sqsubseteq g_{1}\sqsubseteq g_{2}\sqsubseteq\ldots$ be an $\omega$-chain in $\mathcal{D}$. Then, by the monotonicity of $-^{\ddagger}$, $\{g^{\ddagger}_{n}\}_{n}$ is also an $\omega$-chain in $\mathcal{D}^{\ddagger}$. We should show that

[TABLE]

for all $a$ from $(\{\bot\}\cup\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty))$. Pick such an $a$. If $a=\bot$, the required holds because both sides become $\bot$. If $a\neq\bot$ but $g_{n}(a)=\bot$ for all $n$, then again both sides of the equality are $\bot$, and so the equality holds. If neither of these two cases hold, there exists $n$ such that $g_{n}(a)\neq\bot$ and $g_{n}(a)=g_{m}(a)$ for all $m\geq n$. Let $(s^{\prime},r^{\prime},w^{\prime},p^{\prime})=g_{n}(a)$, and $(s,r)=a$. Then,

[TABLE]

∎

Lemma 4.6. For every command $C$, its semantics ${\llbracket C\rrbracket}_{d}$ is well-defined and belongs to $\mathcal{D}$.

Proof.

The proof that ${\llbracket C\rrbracket}_{d}$ is well-defined and measurable proceeds by induction on the structure of commands $C$.

• •

Case of ${\llbracket\mathsf{{skip}}\rrbracket}_{d}$:

The semantics is clearly well-defined in this case. Note that ${\llbracket\mathsf{{skip}}\rrbracket}_{d}$ is essentially the pairing of the identity function and a constant function. Since both of these functions are measurable and the pairing preserves measurability, ${\llbracket\mathsf{{skip}}\rrbracket}_{d}$ is measurable. Also, since ${\llbracket\mathsf{{skip}}\rrbracket}_{d}$ neither changes nor depends on the $r$ component, it is local.

• •

Case of ${\llbracket x:=E\rrbracket}_{d}$:

The proof is identical to the previous case except that instead of the identity function, we should use the function $(s,r)\longmapsto(s[x\mapsto{\llbracket E\rrbracket}s],r)$, which is measurable because it can be constructed by composing and pairing three measurable functions, ${\llbracket E\rrbracket}$, $(s,v)\longmapsto s[x\mapsto v]$, and projection.

• •

Case of ${\llbracket\mathsf{{if}}\,B\,\{C_{0}\}\,\mathsf{{else}}\,\{C_{1}\}\rrbracket}_{d}$:

The well-definedness and the locality are immediate consequences of the induction hypothesis on ${\llbracket C_{0}\rrbracket}_{d}$ and ${\llbracket C_{1}\rrbracket}_{d}$. Let $A\subseteq\{\bot\}\cup\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty)$ be a measurable set. Then,

[TABLE]

Moreover, $\{\mathit{{true}}\}$ is measurable, and so is its complement. Since ${\llbracket B\rrbracket}$ is a measurable function, $\{(s,r)\mid{\llbracket B\rrbracket}s=\mathit{{true}}\}$ and $\{(s,r)\mid{\llbracket B\rrbracket}s\not=\mathit{{true}}\}$ are measurable. As a result the above set is measurable. So, ${\llbracket\mathsf{{if}}\,B\,\{C_{0}\}\,\mathsf{{else}}\,\{C_{1}\}\rrbracket}_{d}$ is well-defined, local and measurable as claimed.

• •

Case of ${\llbracket C_{0};C_{1}\rrbracket}_{d}$:

By induction hypothesis, both ${\llbracket C_{0}\rrbracket}_{d}$ and ${\llbracket C_{1}\rrbracket}_{d}$ are well-defined, local and measurable. This and Lemma A.3 imply that ${\llbracket C_{1}\rrbracket}^{\ddagger}_{d}$ is also local and measurable. The function ${\llbracket C_{0};C_{1}\rrbracket}_{d}$ is nothing but the composition of two well-defined measurable functions ${\llbracket C_{1}\rrbracket}^{\ddagger}_{d}$ and ${\llbracket C_{0}\rrbracket}_{d}$, and so it satisfies the claimed well-definedness and measurability. To see the locality of ${\llbracket C_{0};C_{1}\rrbracket}_{d}$, suppose that

[TABLE]

Then, there exist $s^{\prime},r^{\prime},w^{\prime},p^{\prime}$ such that

[TABLE]

By the locality of ${\llbracket C_{0}\rrbracket}_{d}$ and ${\llbracket C_{1}\rrbracket}^{\ddagger}_{d}$, there exist $r^{\prime}_{0}$ and $r^{\prime\prime}_{0}$ such that

[TABLE]

Note that $r_{0}^{\prime}\#r_{0}^{\prime\prime}$. Let $r_{0}\triangleq r_{0}^{\prime}\uplus r_{0}^{\prime\prime}$. Then,

[TABLE]

The first equality from above uses the locality of ${\llbracket C_{0}\rrbracket}_{d}$. Thus, the first conjunct in the definition of locality holds if we use $r_{0}$ as a witness in the conjunct. To prove the second in the definition, consider $r_{1}$ such that $r_{1}\#r$. Then, $r_{1}\#r^{\prime}$. Thus, by the locality of ${\llbracket C_{0}\rrbracket}_{d}$ and ${\llbracket C_{1}\rrbracket}^{\ddagger}_{d}$, we have

[TABLE]

• •

Case of ${\llbracket\mathsf{{while}}\,B\,\{C\}\rrbracket}_{d}$:

Since $\mathcal{D}$ is an $\omega$-cpo and has the least element, we simply need to prove that $G$ in the definition of ${\llbracket\mathsf{{while}}\,B\,\{C\}\rrbracket}_{d}$ is well-defined and continuous so that Kleene’s fixed point theorem applies (Lemma 4.5). By induction hypothesis, ${\llbracket C\rrbracket}_{d}$ is well-defined, local, and measurable. Thus, $G$ is a well-defined function from $\mathcal{D}$ to $\mathcal{D}$; the proof that the image of $G$ consists of local measurable functions is similar to the case of the if statement extended with the argument with the case of sequencing. Furthermore, $G$ is monotone because both the lifting $-^{\ddagger}$ and the function composition are monotone. Thus, $G$ transforms an $\omega$-chain into an $\omega$-chain. Finally, $G$ preserves the least upper bound of such an $\omega$-chain. This is because both the lifting $-^{\ddagger}$ and the function composition preserve the least upper bound. We have just shown that $G$ is well-defined and continuous.

• •

Case of ${\llbracket x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2})\rrbracket}_{d}$:

The well-definedness and the locality are immediate. To prove the measurability, pick a measurable subset $A$ of $\{\bot\}\cup(\mathit{{Store}}\times\mathit{{RDB}}\times[0,\infty)\times[0,\infty)$. Then,

[TABLE]

where

[TABLE]

We show that all of $A_{0}$, $A_{1}$, $A_{2}$ and $A_{3}$ are measurable. The measurability of the first three sets follow from the measurability of the following two sets:

[TABLE]

The set $A^{\prime}_{0}$ is measurable because it is the inverse image of the composition of two measurable functions (namely, ${\llbracket E_{1}\rrbracket}$ and projection) to the measurable set $(0,\infty)$. For the measurability of $A^{\prime}_{1}$, note that

[TABLE]

The set on the RHS is measurable, because it is a countable union of the intersection of two measurable subsets; the measurability of $\{(s,r)\,\mid\,\mathrm{dom}(r)=K\}$ comes from the fact that $\mathit{{RDB}}$ is built by the disjoint-union construction over $[K\to\mathbb{R}]$ for all finite $K$, and the measurability of $\{(s,r)\,\mid\,{\llbracket S\rrbracket}s\in K\}$ follows from the measurability of the function ${\llbracket S\rrbracket}$ and the set $K$. It remains to show the measurability of $A_{3}$. To do so, we note that $A_{3}$ is the inverse image of a function constructed by applying measurability-preserving operators on ${\llbracket S\rrbracket}$, ${\llbracket E_{1}\rrbracket}$, ${\llbracket E_{2}\rrbracket}$ and $\mathcal{N}$, and the application and restriction functions on random databases. Since the first four are measurable, it is sufficient to show the measurability of the last two, which we spell out below:

[TABLE]

For any measurable subset $A$ of $\mathbb{R}$,

[TABLE]

All the explicitly-defined sets on the RHS of this equation are measurable, and the RHS only uses finite or countable union. Thus, $f_{1}^{-1}(A)$ is measurable. Now consider a measurable subset $A^{\prime}$ of $\mathit{{RDB}}$. Then,

[TABLE]

Again all the explicitly-defined sets on the RHS of this equation are measurable, and the RHS only uses finite or countable union. Thus, $f_{2}^{-1}(A^{\prime})$ is measurable.

• •

Case of ${\llbracket\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})\rrbracket}_{d}$:

The well-definedness and the locality are immediate in this case. Pick a measurable set

[TABLE]

Then,

[TABLE]

Both sets in the union here are measurable because all of $\mathcal{N}$, ${\llbracket E_{0}\rrbracket}$, ${\llbracket E_{1}\rrbracket}$ and ${\llbracket E_{2}\rrbracket}$ are measurable and the composition and the pairing operators on functions preserve measurability. Thus, the union itself is measurable.

∎

Lemma 4.7. For all $g\in\mathcal{D}$, the following functions from $\mathit{{Store}}\times\mathit{{RDB}}$ to $\mathbb{R}$ and $\{\bot\}\cup\mathit{{Store}}$ are measurable: $(s,r)\longmapsto\mathit{{dens}}(g,s)(r)$ and $(s,r)\longmapsto\mathit{{get}}(g,s)(r)$.

Proof.

We have to show that for every real $v\in\mathbb{R}$ and measurable subset $A\subseteq\{\bot\}\cup\mathit{{Store}}$,

[TABLE]

are measurable. The set $B_{0}$ is measurable because

[TABLE]

and the three sets on the RHS of the equation are all measurable. The measurability of the other $B_{1}$ follows from a similar argument based on case split. That is,

[TABLE]

and all three sets on the RHS of the equation are measurable. Thus, $B_{1}$ is a measurable subset. ∎

In the next lemma, we regard $\{\bot\}\cup\mathit{{Store}}$ as a partially-ordered set with the expected order: $a\sqsubseteq b\iff a=\bot\vee a=b$, for all $a,b\in\{\bot\}\cup\mathit{{Store}}$.

Lemma A.4.

Let $\{g_{n}\}_{n}$ be an $\omega$-chain in $\mathcal{D}$. Then, for all stores $s$ and random databases $r$, the sequences $\{\mathit{{dens}}(g_{n},s)(r)\}_{n}$ and $\{\mathit{{get}}(g_{n},s)(r)\}_{n}$ are increasing,

[TABLE]

Proof.

Let $s$ be a store, $r$ be a random database, and $\{g_{n}\}_{n}$ be an $\omega$-chain in $\mathcal{D}$.

We first show that the sequences $\{\mathit{{dens}}(g_{n},s)(r)\}_{n}$ and $\{\mathit{{get}}(g_{n},s)(r)\}_{n}$ are increasing. We do so by showing the monotonicity of $\mathit{{dens}}(-,s)(r)$and $\mathit{{get}}(-,s)(r)$. Consider $g,g^{\prime}\in\mathcal{D}$ such that $g\sqsubseteq g^{\prime}$. Then, $g(s,r)=\bot$ or $g(s,r)=g^{\prime}(s,r)$. In the former case, $\mathit{{dens}}(g,s)(r)=0\leq\mathit{{dens}}(g^{\prime},s)(r)$ and $\mathit{{get}}(g,s)(r)=\bot\sqsubseteq\mathit{{get}}(g^{\prime},s)(r)$. In the latter case, $\mathit{{dens}}(g,s)(r)=\mathit{{dens}}(g^{\prime},s)(r)$ and $\mathit{{get}}(g,s)(r)=\mathit{{get}}(g^{\prime},s^{\prime})(r)$. Thus, $\mathit{{dens}}(-,r)(s)$ and $\mathit{{get}}(-,r)(s)$ are monotone.

Next we show that $\mathit{{dens}}(-,s)(r)$ and $\mathit{{get}}(-,s)(r)$ preserve the limit. By what we have already shown, the sequence $\{g_{n}(s,r)\}_{n}$ is monotone. But this sequence is ultimately stable. The stable value becomes the least upper bound $\bigsqcup_{n}g_{n}(s,r)$, which is also equal to $(\bigsqcup_{n}g_{n})(s,r)$. Thus, $\{\mathit{{dens}}(g_{n},s)(r)\}_{n}$ is ultimately stable as well, and its limit is

[TABLE]

By the same reason, $\{\mathit{{get}}(g_{n},s)(r)\}_{n}$ is ultimately stable, and has the limit

[TABLE]

∎

Lemma 4.8. For all non-negative bounded measurable functions $h:(\{\bot\}\cup\mathit{{Store}})\times\mathit{{RDB}}\to\mathbb{R}$, stores $s$, and functions $g_{1},g_{2}\in\mathcal{D}$, we have that

[TABLE]

Proof.

For $g\in\mathcal{D}$ and $s\in\mathit{{Store}}$, let

[TABLE]

Here $r\setminus K$ means the restriction of $r$ to $\mathrm{dom}(r)\setminus K$. Both $R_{1}(g,s)$ and $R_{2}(g,s)$ are measurable. The measurability of the latter is an immediate consequence of the measurability of $g$. For the former, we note that for every measurable subset $A$ of $\mathit{{RDB}}$ with $[]\not\in A$,

[TABLE]

The partially applied function $g(s,-)$ is measurable, so that the RHS of this equation is the countable union of measurable subsets. If $A=\{[]\}$, then

[TABLE]

Thus, $R_{1}(g,s)^{-1}(A)$ is measurable. The remaining case is that $[]\in A$. We can be handle this case by splitting $A$ into $\{[]\}$ and $A\setminus\{[]\}$, and dealing with the slit cases separately.

For $g\in\mathcal{D}$, $s\in\mathit{{Store}}$, and $r\in\mathit{{RDB}}$, let $\mathit{{consumed}}(g,s,r)$ be the predicate defined by

[TABLE]

Then, $r\longmapsto{\mathds{1}_{[{\mathit{{consumed}}(g,s,r)}]}}$ is a measurable function.

Now we note a few useful equalities about the entities $h,s,g_{1},g_{2}$ assumed in the lemma.

First, for all $r$, if

[TABLE]

then

[TABLE]

The first equality is just an immediate consequence of the definitions of $R_{1}$ and $R_{2}$, and the other two equalities hold mainly because $g_{1}$ and $g_{2}$ are local.

Next, consider the density functions $f_{1}$ on $\mathit{{RDB}}$ and $f_{2}$ on $\mathit{{RDB}}\times\mathit{{RDB}}$:

[TABLE]

The densities are taken with respect to $\rho$ and $\rho\otimes\rho$, respectively. There is a measurable bijection $\beta$ between the supports of $f_{1}$ and $f_{2}$, and the pushforward of the measure of $f_{1}$ by this bijection is the measure of $f_{2}$. Concretely, the function $\beta$ is:

[TABLE]

where $r_{d}\in\mathit{{RDB}}$ is $[``s"\mapsto\mathit{{nil}}]$. Note that $\beta$ is not a bijection. But the restriction of $\beta$ to $A_{1}=\{r\,\mid\,f_{1}(r)>0\}$ is injective, and its image is $A_{2}=\{(r_{1},r_{2})\,\mid\,f_{2}(r_{1},r_{2})>0\}$. These properties follow from the definitions of $\mathit{{consumed}}$, $R_{1}$ and $R_{2}$ and the locality of $g$. The function $\beta$ is measurable because $R_{1}(g_{1},s)$ and $R_{2}(g_{2},s)$ are measurable functions and $\{r\,\mid\,\mathit{{consumed}}(g_{2}^{\ddagger}\circ g_{1},s,r)\}$ is a measurable set. It remains to show the measure preservation of $\beta$. Pick a measurable subset $A^{\prime}$ of $A_{2}$. We have to show that

[TABLE]

We calculate the equation as follows:

[TABLE]

Using what we have shown so far, we calculate the claim of the lemma. In the calculation we often omit the $g_{1}$ and $s$ parameters from $R_{1}(g_{1},s)(r)$ and $R_{2}(g_{1},s)(r)$, and just write $R_{1}(r)$ and $R_{2}(r)$.

[TABLE]

∎

Theorem 4.9. For all non-negative bounded measurable monotone functions $h:(\{\bot\}\cup\mathit{{Store}})\times\mathit{{RDB}}\to\mathbb{R}$ and states $\sigma$,

[TABLE]

Proof.

We prove the theorem by induction on the structure of $C$.

When $C\equiv\mathsf{{skip}}$, we derive the claimed equality as follows:

[TABLE]

Next we handle the case that $C\equiv(x:=E)$.

[TABLE]

We move on to the case that $C\equiv(\mathsf{{if}}\,B\,\{C_{0}\}\,\mathsf{{else}}\,\{C_{1}\})$, and prove the desired equality:

[TABLE]

Now we prove the claimed equality for $C\equiv(C_{0};C_{1})$.

[TABLE]

The third and fifth equalities follow from the induction hypothesis, and the last equality from Lemma 4.8.

The next case is $C\equiv(\mathsf{{while}}\,B\,\{C_{0}\})$. Let $\mathcal{R}\subseteq\mathcal{K}\times\mathcal{D}$ be the relation defined by

[TABLE]

Note that the least elements from $\mathcal{K}$ and $\mathcal{D}$ are related by $\mathcal{R}$ because they make both sides of the equality in the definition of $\mathcal{R}$ be zero. Furthermore, for all $\omega$-chains $\{\kappa_{n}\}_{n}$ and $\{g_{n}\}_{n}$ in $\mathcal{K}$ and $\mathcal{D}$, if $\kappa_{n}\,[\mathcal{R}]\,g_{n}$ for all $n$, then

[TABLE]

This is because for all stores $\sigma$ and non-negative bounded measurable functions $h$,

[TABLE]

The first equation follows from Lemma A.1, and the third from the monotone convergence theorem and Lemma A.4. The last equality holds because of Lemma A.4 again.

Thus, we can complete the proof of this case if we show one more thing. Let $K$ and $G$ be the functions on $\mathcal{K}$ and $\mathcal{D}$ used in the interpretation of the loop (namely, $\mathsf{{while}}\,B\,\{C_{0}\}$) in the two semantics. What we have to show is:

[TABLE]

Pick $\kappa$ and $g$ with $\kappa\,[\mathcal{R}]\,g$. Then, for all $\sigma$ and non-negative bounded measurable $h$,

[TABLE]

But

[TABLE]

The last equality uses Lemma 4.8. Thus, if we continue our calculation that we paused momentarily, we get

[TABLE]

Next we handle the case that $C\equiv(x:=\mathsf{{sample}}_{\mathsf{{norm}}}(S,E_{1},E_{2}))$.

[TABLE]

It remains to prove the case that $C\equiv\mathsf{{score}}_{\mathsf{{norm}}}(E_{0},E_{1},E_{2})$. Here is our proof for this last case:

[TABLE]

∎

Appendix B Proof of Theorem in §5

Theorem 5.1. Let $C$ be a model, $D_{\theta}$ be a guide, and $N\neq 0\in\mathbb{N}$. Define $\mathrm{KL}_{(-)}:\mathbb{R}^{p}\to\mathbb{R}_{\geq 0}$ as $\mathrm{KL}_{\theta}$ $\triangleq$ $\mathrm{KL}(\mathit{{dens}}(D_{\theta},s_{I}){||}\mathit{{dens}}(C,s_{I})/Z_{C})$. Then, $\mathrm{KL}_{(-)}$ is well-defined and continuously differentiable with

[TABLE]

if

1. (R1)

$\mathit{{dens}}(C,s_{I})(r)=0\implies\mathit{{dens}}(D_{\theta},s_{I})(r)=0$, for all $r\in\mathit{{RDB}}$ and $\theta\in\mathbb{R}^{p}$; 2. (R2)

for all $(r,\theta,j)\in\mathit{{RDB}}\times\mathbb{R}^{p}\times[p]$, the function $v\longmapsto\mathit{{dens}}(D_{\theta[j:v]},s_{I})(r)$ on $\mathbb{R}$ is differentiable; 3. (R3)

for all $\theta\in\mathbb{R}^{p}$,

[TABLE] 4. (R4)

for all $(\theta,j)\in\mathbb{R}^{p}\times[p]$, the function

[TABLE]

on $\mathbb{R}$ is continuously differentiable; 5. (R5)

for all $\theta\in\mathbb{R}^{p}$,

[TABLE] 6. (R6)

for all $\theta\in\mathbb{R}^{p}$,

[TABLE]

As we explained in the main text already, the theorem (except the continuous differentiability of $\mathrm{KL}_{\theta}$) is well-known with a well-known proof. We include its proof for completeness and also to help the reader to see why each of the requirements in the theorem is needed.

Proof.

$\mathrm{KL}_{\theta}$ is well-defined because of the following derivation with R1 and R3:

[TABLE]

It remains to show the continuous differentiability of $\mathrm{KL}_{\theta}$ and Equation (6). Recall the fact that for any function $f:\mathbb{R}^{N}\to\mathbb{R}^{M}$, if all partial derivatives of $f$ exist and are continuous, $f$ is continuously differentiable. Because of the fact and the linearity of expectation, it suffices to show that all partial derivatives of $\mathrm{KL}_{\theta}$ are well-defined and continuous, and

[TABLE]

Here is a well-known derivation of (20):

[TABLE]

The well-definedness and the continuity of $(\nabla_{\theta}\mathrm{KL}_{\theta})_{j}$ for all $j\in[p]$ follow from (21) and R4, which concludes the proof.

Note that R2 is used to guarantee the well-definedness of $\nabla_{\theta}\log\mathit{{dens}}(D_{\theta},s_{I})(r)$ in (20); R5 and R6 are used in (22) and (23), respectively. ∎

Appendix C Proofs of Examples, Lemma, and Theorems in §6

Theorem 6.1. Under our assumption in §6.1,

[TABLE]

Proof.

For $i\in[M]$ and $r\in A_{i}$, define

[TABLE]

Then,

[TABLE]

The last inequality uses a well-known result that the differential entropy of any multivariate normal distribution is finite. ∎

Example 6.2. Consider guides $D_{(i,\theta)}$ defined as follows ($i=1,2$):

[TABLE]

where for some $n\geq 1$ and $c\neq 0\in\mathbb{R}$,

[TABLE]

Then, the entropies of $\mathit{{dens}}(D_{(i,\theta)},s_{I})$’s are all undefined.

Proof.

We prove that

[TABLE]

by bounding the the expectation as follows.

[TABLE]

Here the two terms in the last equation are bounded as follows.

[TABLE]

Note that

[TABLE]

because $\int_{[0,\epsilon]}1/v^{n}\,\mathrm{d}v=\infty$ for any $\epsilon>0$ and $n\geq 1$, and $\mathbb{E}_{\mathcal{N}(v;\mu,\sigma)}[\exp(c\cdot v^{3})]=\infty$ for any $\mu\in\mathbb{R}$ and $\sigma\in\mathbb{R}_{>0}$. Hence, we obtain the desired result. ∎

Theorem 6.3. Under our assumption in §6.1, the condition (13) implies

[TABLE]

and thus the requirement R3 (i.e., the objective in (5) is well-defined).

Proof.

For $i\in[M]$, $j\in[N_{i}]$, $\alpha\in K_{i}$, and $r\in A_{i}$, let

[TABLE]

We note that

[TABLE]

Thus, we can have

[TABLE]

if the following two inequalities hold:

[TABLE]

We show that the inequalities are indeed true. By the condition (13) and the parts (a) and (b) of Lemma 6.6 (full version), there exist functions $f_{i},u_{i}\in\mathcal{A}_{i}$ such that for every $r\in A_{i}$,

[TABLE]

For every $f\in\mathcal{A}_{i}$, by the part (e) of Lemma 6.6 (full version), we have that

[TABLE]

The inequalities in (25) follow from the properties in (26) and (27).

By an argument similar to what we have just given, we also have that

[TABLE]

Using (24) and (28), we calculate the claim of the theorem:

[TABLE]

The last inequality follows from (24) and (28). Finally, the requirement R3 follows from the above and Theorem 6.1:

[TABLE]

∎

Example 6.4. Consider models $C_{1},\ldots,C_{4}$ and a guide $D_{\theta}$ defined as follows:

[TABLE]

where for some $n\geq 1$ and $c\neq 0\in\mathbb{R}$,

[TABLE]

Then, the requirement R3 does not hold (i.e., the objective in (5) is undefined) for all $i=1,\ldots,4$.

Proof.

Since we have

[TABLE]

it suffices to show that

[TABLE]

Case of $i=1,2$: The quantity in (29) is bounded as follows.

[TABLE]

Here the two terms in the last equation are bounded as follows.

[TABLE]

We bound the last term in the last equation as follows, assuming $c>0$.

[TABLE]

Here the last inequality comes from that the two expectations over $x_{2}$ are finite while the other two over $x_{1}$ are both $\infty$, because $\int_{[0,\epsilon]}1/v^{n}\,\mathrm{d}v=\infty$ for any $\epsilon>0$ and $n\geq 1$, and $\mathbb{E}_{\mathcal{N}(v;\mu,\sigma)}[\exp(c\cdot v^{3})]=\infty$ for any $\mu\in\mathbb{R}$ and $\sigma\in\mathbb{R}_{>0}$. Note that we can use a similar calculation for the case of $c<0$. This proves the equality in (29).

Case of $i=3,4$: The proof is similar to the above case, so we omit it. Remark that in the proof, we additionally use the following fact: $\int_{[0,1]}|\log v|\,\mathrm{d}v<\infty$. ∎

Lemma 6.6 (Full Version). Pick $i\in[M]$. Let $(\alpha_{1},\ldots,\alpha_{J})$ be an enumeration of the elements in $K_{i}$, and $\overline{r}\triangleq(r(\alpha_{1}),\ldots,r(\alpha_{J}))\in\mathbb{R}^{J}$ for $r\in A_{i}$.

• (a)

For every $l_{1},l_{2}\in\mathcal{A}_{i}$ and every $\odot\in\{+,-,\times,/,\max\}$, there exists $l\in\mathcal{A}_{i}$ such that

[TABLE]

• (b)

For every polynomial $\mathit{poly}:\mathbb{R}^{J}\to\mathbb{R}$, there exists $l\in\mathcal{A}_{i}$ such that

[TABLE]

In particular, for every $l\in\mathcal{A}_{i}$, there exists $l^{\prime}\in\mathcal{A}_{i}$ such that

[TABLE]

• (c)

For every affine-bounded neural network $\mathit{{nn}}:\mathbb{R}^{J}\to\mathbb{R}$, there exists $l\in\mathcal{A}_{i}$ such that

[TABLE]

• (d)

There exist affine functions $l_{1},l_{2}:\mathbb{R}\to\mathbb{R}$ such that

[TABLE]

• (e)

For every $l\in\mathcal{A}_{i}$ and every $J$-dimensional normal distribution $q$,

[TABLE]

Proof of (a).

We split the cases.

• •

Case $\odot\in\{\times,/\}$: Take $l=l_{1}\pm l_{2}$.

• •

Case $\odot=\max$: For $m=1,2$, let

[TABLE]

where $c_{m,j}\in\mathbb{R}$. For all $r\in A_{i}$, we have

[TABLE]

where

[TABLE]

Hence, for every $r\in A_{i}$,

[TABLE]

• •

Case $\odot\in\{+,-\}$: Using the case $\odot=\max$, we have the following for some $l\in\mathcal{A}_{i}$:

[TABLE]

Proof of (b). Let $n$ be an index in $\{1,\ldots,J\}$ such that

[TABLE]

for some polynomial $\mathit{poly}^{\prime}:\mathbb{R}^{n}\to\mathbb{R}$. We show that there exists $l\in\mathcal{A}_{i}$ satisfying

[TABLE]

Note that the claim of the lemma follows from this. Our proof is by induction on $n$.

• •

Case $n=1$: Let $v\in\mathbb{R}$ and $\mathit{poly}^{\prime}(v)=\sum_{d=0}^{D}c_{d}\cdot v^{d}$ where $c_{d}\in\mathbb{R}$. We use the following fact.

Fact: For every $d\in\mathbb{N}$, there exists an affine function $l_{d}:\mathbb{R}\to\mathbb{R}$ such that

[TABLE]

Proof.

Since $v^{d}=\mathcal{O}(\exp(v))$ for positive $v\in\mathbb{R}$, there exist $C_{1},C_{2}>0$ such that $|v|>C_{1}$ implies $|v|^{d}<C_{2}\exp(|v|)$. Since $|v|^{d}$ is continuous, $C_{3}\triangleq\sup\{|v|^{d}\,\mid\,|v|\leq C_{1}\}$ satisfies $0\leq C_{3}<\infty$. Putting these together, for every $v\in\mathbb{R}$,

[TABLE]

where $C_{4}\triangleq\max(C_{2},1)\cdot(C_{3}+1)\geq 1$. Thus, $l_{d}(v)\triangleq(v+\log C_{4})$ is a desired affine function. ∎

Using the fact, we have

[TABLE]

for some $l^{\prime}_{d},l\in\mathcal{A}_{i}$. The existence of $l$ satisfying the last inequality follows from (a).

• •

Case $n>1$: Let $v\in\mathbb{R}^{n}$, and define $v_{1:n-1}$ to be the projection of $v$ to its first $n-1$ components. Also, let

[TABLE]

where $\mathit{poly}^{\prime}_{d}:\mathbb{R}^{n-1}\to\mathbb{R}$ is a polynomial. Using the induction hypothesis, we have

[TABLE]

for some affine functions $l_{d}:\mathbb{R}\to\mathbb{R}$, and functions $l^{\prime}_{d},l^{\prime\prime}_{d},l\in\mathcal{A}_{i}$. The last inequality follows from (a).

Proof of (c). By the definition of affine-bounded neural network, there exist functions

[TABLE]

and affine functions

[TABLE]

such that (i) $n_{1}=J$ and $n_{d+1}=1$; (ii) $\mathit{{nn}}=f_{d}\circ\cdots\circ f_{1}$; and (iii) for all $1\leq j\leq d$ and $v\in\mathbb{R}^{n_{j}}$,

[TABLE]

We prove the claim of the lemma by induction on the depth $d$. For $d=1$, the claim follows immediately from the defining property (iii) of the affine-bounded neural network $\mathit{nn}$. Assume $d>1$. By induction hypothesis, there exist affine functions $l_{1},\ldots,l_{n_{d}}:\mathbb{R}^{J}\to\mathbb{R}$ such that for all $k\in\{1,\ldots,n_{d}\}$ and $v\in\mathbb{R}^{J}$,

[TABLE]

This is because the $(f_{d-1}\circ\cdots\circ f_{1})(v)_{k}$ are affine-bounded neural networks of depth $d-1$. Let

[TABLE]

Using what we have proved so far, we calculate:

[TABLE]

Let $l(r)=c_{0}+\sum_{k=1}^{n_{d}}|c_{k}|\cdot l_{k}(|r(\alpha_{1})|,\ldots,|r(\alpha_{J})|)$. The function $l$ is essentially the composition of two affine functions, and so it belongs to $\mathcal{A}_{i}$. This and the conclusion of the above calculation show that $l$ is the function in $\mathcal{A}_{i}$ that we look for.

Proof of (d). Since $\log(1+v)\geq v/2$ for all $v\in[0,1]$, we have

[TABLE]

Let $l_{1}(v)=-v-\log 2$ and $l_{2}(v)=v+\log 2$. As $\mathrm{softplus}(-)$ is increasing, we have

[TABLE]

Proof of (e). Let $l(r)=c_{0}+\sum_{j=1}^{J}c_{j}\cdot|r(\alpha_{j})|$, where $c_{j}\in\mathbb{R}$. The moment generating function $G_{p}:\mathbb{R}^{J}\to\mathbb{R}$ of a random variable of density $p(v)$ on $\mathbb{R}^{J}$ is defined by $G_{p}(t)\triangleq\mathbb{E}_{p(v)}\left[{\exp(t\cdot v)}\right]$ where $t\cdot v$ is the inner product of $t$ and $v$. It is well-known that the moment generating functions are well defined for multivariate normal distributions $q(v)$ on $\mathbb{R}^{J}$ in the lemma, i.e., $G_{q}(t)<\infty$ for all $t\in\mathbb{R}^{J}$. From this fact and that $\exp(|v|)\leq\exp(v)+\exp(-v)$, it follows that

[TABLE]

∎

The proof of Theorem 6.7 uses two more ingredients: Theorem C.1 which, compared to Theorem 6.8, additionally guarantees the continuity of $\nabla_{v}\!\int\!f_{v}\,\mathrm{d}\mu$ by assuming the continuity of $\nabla_{v}f_{v}$; and Lemma C.2 which bounds the densities of multiple normal distributions using the density of a single, fixed normal distribution.

Theorem C.1.

Let $V\subset\mathbb{R}$ be an open interval, and $(X,\Sigma,\mu)$ be a measure space. Suppose that a measurable function $f:V\times X\to\mathbb{R}$ satisfies the following conditions:

• (a)

The function $v\in V\longmapsto\int\mu(\mathrm{d}x)\,f_{v}(x)\in\mathbb{R}$ is well-defined.

• (b)

For almost all $x\in X$ (w.r.t. $\mu$), the function $v\in V\longmapsto\nabla_{v}f_{v}(x)\in\mathbb{R}$ is well-defined, continuous.

• (c)

There is a measurable function $h:X\to\mathbb{R}$ such that $\int\mu(\mathrm{d}x)\,h(x)$ is well-defined and $\left|\nabla_{v}f_{v}(x)\right|\leq h(x)$ for all $v\in V$ and almost all $x\in X$ (w.r.t. $\mu$).

Then, for all $v\in V$, both sides of the below equation are well-defined, and the equality holds:

[TABLE]

Moreover, the function $v\in V\longmapsto\nabla_{v}\int\mu(\mathrm{d}x)\,f_{v}(x)$ is continuous.

Proof.

Theorem 6.8 implies the former part of the conclusion: for all $v\in V$, each side of (30) is well-defined and the equation (30) holds. To prove the latter part of the conclusion, define $G:V\to\mathbb{R}$ and $g:V\times X\to\mathbb{R}$ by

[TABLE]

and let $(v_{n})_{n\in\mathbb{N}}$ be any sequence in $V$ such that $\lim_{n\to\infty}v_{n}=v^{*}\in V$. Then, it suffices to show that

[TABLE]

To apply the dominated convergence theorem, define $g_{n},\,g^{*}:X\to\mathbb{R}$ by $g_{n}(x)\triangleq g(v_{n},x)$ and $g^{*}(x)\triangleq g(v^{*},x)$. By (b), the function $g(\cdot,x):V\to\mathbb{R}$ is continuous for almost all $x\in X$ (w.r.t. $\mu$), so

[TABLE]

By (c), for all $n\in\mathbb{N}$, $|g_{n}(x)|\leq h(x)$ for almost all $x\in X$ (w.r.t. $\mu$), and $h(x)$ is integrable (w.r.t. $\mu$). Therefore, we can apply the dominated convergence theorem to $(g_{n})_{n\in\mathbb{N}}$ and $g^{*}$ with $h$, which gives:

[TABLE]

∎

Lemma C.2.

Suppose $[\mu_{1},\mu_{2}]\subset\mathbb{R}$ and $[\sigma_{1},\sigma_{2}]\subset(0,\infty)$ Then, there exists $c\in\mathbb{R}$ such that

[TABLE]

for all $v\in\mathbb{R}$, $\mu\in[\mu_{1},\mu_{2}]$, and $\sigma\in[\sigma_{1},\sigma_{2}]$.

Proof.

We claim that $c=\frac{2\sigma_{2}}{\sigma_{1}}\exp(\frac{1}{6\sigma_{1}^{2}}\max(\mu_{1}^{2},\mu_{2}^{2}))$ is the desired coefficient. For the constant $c$, the claim we need show is that

[TABLE]

for any $v\in\mathbb{R}$, $\mu\in[\mu_{1},\mu_{2}]$, and $\sigma\in[\sigma_{1},\sigma_{2}]$. Since $\frac{1}{\sqrt{2\pi}\sigma_{1}}\geq\frac{1}{\sqrt{2\pi}\sigma}$, it suffices to show

[TABLE]

The last inequality holds because

[TABLE]

∎

Theorem 6.7. If both our assumption in §6.1 and the condition (13) hold, then the condition (17) implies the requirements R4-R6.

Proof.

Before staring the proof, we note that by the condition (17), there is an open interval $V_{(i,\alpha,j,\theta)}\subseteq\mathbb{R}$ for all $i\in[M]$, $\alpha\in K_{i}$, $j\in[p]$ and $\theta\in\mathbb{R}^{p}$, such that

• •

$V_{(i,\alpha,j,\theta)}$ contains $\theta_{j}$;

• •

the function $v\in\mathbb{R}\longmapsto\mu_{(i,\alpha)}(\theta[j:v])$ is continuously differentiable in $V_{(i,\alpha,j,\theta)}$;

• •

the function $v\in\mathbb{R}\longmapsto\sigma_{(i,\alpha)}(\theta[j:v])$ is continuously differentiable in $V_{(i,\alpha,j,\theta)}$, and $\inf\{\sigma_{(i,\alpha)}(\theta[j:v])\,|\,v\in V_{(i,\alpha,j,\theta)}\}>0$.

Let $F_{1}$ and $F_{2}$ be functions from $\mathbb{R}^{p}\times\mathit{{RDB}}$ to $\mathbb{R}$ defined as follows: for all $\theta\in\mathbb{R}^{p}$ and $r\in\mathit{{RDB}}$,

[TABLE]

We have to show that for every $\theta\in\mathbb{R}^{p}$ and $j\in[p]$,

• •

the partial derivatives

[TABLE]

are well-defined and are continuous functions of $\theta_{j}\in\mathbb{R}$;

• •

these derivatives commute with integration:

[TABLE]

We discharge these proof obligations using Theorem C.1. Pick $\theta^{*}\in\mathbb{R}^{p}$ and $j\in[p]$. We instantiate Theorem C.1 twice for $F_{1}$ and $F_{2}$ as follows.

• •

For both cases of $F_{1}$ and $F_{2}$, we set $V$ to be an open interval $(a,b)$ such that both $a$ and $b$ are finite, the closure $[a,b]$ is contained in $V_{(j,\theta^{*})}\triangleq\bigcap_{i\in[M]}\bigcap_{\alpha\in K_{i}}V_{(i,\alpha,j,\theta^{*})}$, and $(a,b)$ contains $\theta^{*}_{j}$. Note that such an open interval exists because $V_{(j,\theta^{*})}$ is open and contains $\theta^{*}_{j}$.

• •

The measurable space $(X,\Sigma)$ in the theorem is $\mathit{{RDB}}$ with its $\sigma$-algebra in both cases.

• •

The measure $\mu$ is the reference measure $\rho$ on $\mathit{{RDB}}$ in both cases.

• •

The measurable function $f$ is instantiated twice for $(v,r)\longmapsto F_{i}(\theta^{*}[j:v],r)$ with $i=1,2$.

Under this instantiation, the conclusion of Theorem C.1 implies all of the requirements R4-R6.

It remains to show that both instantiations satisfy the three provisos of Theorem C.1 marked (a), (b) and (c). We tackle the provisos one-by-one. We use the notations in §6.1 without recalling them explicitly. For $i\in[M]$, $l\in[N_{i}]$, $\alpha\in K_{i}$, $v\in\mathbb{R}$, $r\in A_{i}$ and $\theta\in\mathbb{R}^{p}$, let

[TABLE]

Proof of (a). Since we assume the condition 13, Theorem 6.3 ensures that the integral $\int\mu(\mathrm{d}r)\,F_{1}(\theta^{*}[j:v],r)$ is well-defined for all $v\in V$. The other integral involving $F_{2}$ is the integral of the density of $D_{\theta^{*}[j:v]}$. Since $D_{\theta^{*}[j:v]}$ is a guide, the integral is well-defined and has the value $1$.

Proof of (b). Note that

[TABLE]

Recall that every $r\in A_{i}$ is a map from $K_{i}$ to $\mathbb{R}$. The normal density $\mathcal{N}(v^{\prime};\mu^{\prime},\sigma^{\prime})$ and its log are smooth functions on $(v^{\prime},\mu^{\prime},\sigma^{\prime})\in\mathbb{R}\times\mathbb{R}\times(0,\infty)$. Also, by the condition (17), the functions

[TABLE]

are continuously differentiable for all $i\in[M]$ and $\alpha\in K_{i}$. Thus, they are also continuously differentiable when we restrict their domains to $V$. One more thing to notice is that neither ${\mathds{1}_{[{r\in A_{i}}]}}$ nor $-\log(h_{i}(r)\cdot k_{i}(r))$ depends on $v$, so that they can be viewed as constant functions on $v$, which are obviously smooth. Now note that both $v\in V\longmapsto F_{1}(\theta^{*}[j:v],r)$ and $v\in V\longmapsto F_{2}(\theta^{*}[j:v],r)$ are obtained from these smooth or continuously differentiable functions by function composition, addition and multiplication. Thus, they are continuously differentiable.

Proof of (c). For all $i\in[M]$, we will construct measurable function $H_{i},H^{\prime}_{i}:A_{i}\to[0,\infty)$ such that

• •

both $\int\rho(\mathrm{d}r)\,({\mathds{1}_{[{r\in A_{i}}]}}\cdot H_{i}(r))$ and $\int\rho(\mathrm{d}r)\,({\mathds{1}_{[{r\in A_{i}}]}}\cdot H^{\prime}_{i}(r))$ are well-defined;

• •

for all $r\in A_{i}$ and $v\in V$,

[TABLE]

• •

for all $r\in A_{i}$ and $v\in V$,

[TABLE]

Once we have such $H_{i}$’s, we can use them to define desired measurable functions $H,H^{\prime}:\mathit{{RDB}}\to\mathbb{R}$ as follows:

[TABLE]

Since both $H_{i}$ and $H^{\prime}_{i}$ are integrable for all $i$, so are $H$ and $H^{\prime}$. We prove this only for $H$ below. The case for $H^{\prime}$ is similar.

[TABLE]

Both $H_{i}$ and $H^{\prime}_{i}$ bound the magnitudes of the partial derivatives involving $g_{(i,\theta^{*}[j:v])}$ for every $i$. Using this, we show that $H$ and $H^{\prime}$ meet the condition about bounding the magnitude of a derivative:

[TABLE]

It remains to construct measurable $H_{i},H^{\prime}_{i}:A_{i}\to\mathbb{R}$ with the required properties for every $i\in[M]$. Pick $i\in[M]$. Note that

[TABLE]

We will build $H_{i}$ and $H^{\prime}_{i}$ by finding good bounds for $|\nabla_{v}g_{(i,\theta^{*}[j:v])}(r)|$, $|\log g_{(i,\theta^{*}[j:v])}(r)|$ and $|\log h_{i}(r)\cdot k_{i}(r)|$ and then combining these bounds. We do so in three steps.

First, we build functions that bound $|\nabla_{v}g_{(i,\theta^{*}[j:v])}(r)|$ and $|\log g_{(i,\theta^{*}[j:v])}(r)|$ over $v\in V$. For every $\alpha\in K_{i}$,

[TABLE]

Thus, we have

[TABLE]

for some continuous functions $G_{(1,\alpha,n)}:V_{(j,\theta^{*})}\to\mathbb{R}$ that do not depend on $r$. The continuity of $G_{(1,\alpha,n)}$ follows from our condition (17). Also,

[TABLE]

for some continuous functions $G_{(2,\alpha,n)}:V_{(j,\theta^{*})}\to\mathbb{R}$ that do not depend on $r$. The continuity of $G_{(2,\alpha,n)}$ follows from the condition (17). Since $\overline{V}\subset V_{(j,\theta^{*})}$, the functions $G_{(m,\alpha,n)}$ are all continuous on the closed (bounded) interval $\overline{V}$, where $\overline{V}$ is the closure of $V$. Thus, by the extreme value theorem,

[TABLE]

Since the term $g_{(i,\theta^{*}[j:v])}(r)$ in (31) still depends on $v$, we bound it over $v\in V$ using Lemma C.2 which is instantiated with

[TABLE]

We check that the assumptions of the lemma are satisfied: since $\overline{V}\subset V_{(j,\theta^{*})}$ and the functions

[TABLE]

are continuous by the condition (17), $\mu_{m}$ and $\sigma_{m}$ ($m=1,2$) are all finite by the extreme value theorem; moreover, $\sigma_{1}$ is positive by the choice of $V_{(i,\alpha,j,\theta^{*})}$. Thus, by Lemma C.2, there exists $c\in\mathbb{R}$ such that

[TABLE]

for all $v\in V$, where $c^{\prime}\triangleq c^{|K_{i}|}$ and $g^{\prime}(r)\triangleq\prod_{\alpha\in K_{i}}\mathcal{N}(r(\alpha);0,2\sigma_{2})$.

Combining what we have proved, we obtain the following bounds: for all $v\in V$,

[TABLE]

Second, we bound $|\log h_{i}(r)\cdot k_{i}(r)|$ by $\exp(l_{3}(r))$ for some $l_{2}\in\mathcal{A}_{i}$.

[TABLE]

Finally, we define $H_{i}$ and $H^{\prime}_{i}$ by combining the bounds that we got so far. Here are our definitions for them:

[TABLE]

Then, by the choice of $l_{1}$ and $l_{2}$, for all $r\in A_{i}$ and $v\in V$,

[TABLE]

Also, we have

[TABLE]

and

[TABLE]

We have just shown that $H_{i}$ and $H^{\prime}_{i}$ are the functions that we are looking for. ∎

Appendix D Proof of Theorem in §7

Theorem 7.1 (Soundness). For all commands $C$, we have ${\llbracket C\rrbracket}_{d}\in\gamma({\llbracket C\rrbracket}^{\sharp})$.

Proof.

We prove the theorem by structural induction on $C$. Every case except the one for the loop follows from the assumption made for an abstract element in $\mathcal{T}^{\sharp}$ or an operator on it that corresponds to $C$. The proof of the case of $C\equiv(\mathsf{{while}}\ E\ \{C_{0}\})$ goes as follows. Let $t_{\mathit{fix}}\triangleq(\mathit{{wfix}}\ T)$. Then, by the assumption on the concretisation, $\gamma(t_{\mathit{fix}})$ is an admissible subset of $\mathcal{D}$. Thus, it suffices to prove that when $G$ is the function in the density semantics of ${\llbracket\mathsf{{while}}\ E\ \{C_{0}\})\rrbracket}_{d}$, the image of $G$ on $\gamma(t_{\mathit{fix}})$ is included in $\gamma(t_{\mathit{fix}})$. To do so, pick $g\in\gamma(t_{\mathit{fix}})$. We have to show that $G(g)\in\gamma(t_{\mathit{fix}})$. By the induction hypothesis and the assumption on the abstract composition,

[TABLE]

Also, ${\llbracket\mathsf{{skip}}\rrbracket}_{d}\in\gamma(\mathit{{skip}}^{\sharp})$. Thus,

[TABLE]

The set membership holds because of the assumption on the abstract conditional operator. ∎