On Nesting Monte Carlo Estimators

Introduction

Monte Carlo (MC) methods are used throughout the quantitative sciences. For example, they have become a ubiquitous means of carrying out approximate Bayesian inference (Doucet et al., 2001; Gilks et al., 1995). The convergence of MC estimation has been considered extensively in the literature (Durrett, 2010). However, the implications arising from the nesting of MC schemes, where terms in the integrand depend on the result of separate, nested, MC estimators, is generally less well known. This paper examines the convergence of such nested Monte Carlo (NMC) methods.

Nested expectations occur in wide variety of problems from portfolio risk management (Gordy and Juneja, 2010) to stochastic control (Belomestny et al., 2010). In particular, simulations of agents that reason about other agents often include nested expectations. Tackling such problems requires some form of nested estimation scheme like NMC.

A common class of nested expectations is doubly-intractable inference problems (Murray et al., 2006; Liang, 2010), where the likelihood is only known up to a parameter-dependent normalizing constant. This can occur, for example, when nesting probabilistic programs (Mantadelis and Janssens, 2011; Le et al., 2016). Some problems are even multiply-intractable, such that they require multiple levels of nesting to encode (Stuhlmüller and Goodman, 2014). Our results can be used to show that changes are required to the approaches currently employed by probabilistic programming systems to ensure consistent estimation for such problems (Rainforth, 2017, 2018).

The expected information gain used in Bayesian experimental design (Chaloner and Verdinelli, 1995) requires the calculation of an entropy of a marginal distribution and therefore the expectation of the logarithm of an expectation. By extension, any Kullback-Leibler divergence where one of the terms is a marginal distribution also involves a nested expectation. Hence, our results have important implications for relaxing mean-field assumptions, or using different bounds, in variational inference (Hoffman and Blei, 2015; Naesseth et al., 2017; Maddison et al., 2017) and deep generative models (Burda et al., 2015; Le et al., 2018).

Certain nested estimation problems can be tackled by pseudo-marginal methods (Beaumont, 2003; Andrieu and Roberts, 2009; Andrieu et al., 2010). These consider inference problems where the likelihood is intractable, but can be estimated unbiasedly. From a theoretical perspective, they reformulate the problem in an extended space with auxiliary variables that are used to represent the stochasticity in the likelihood computation, enabling the problem to be expressed as a single expectation.

Our work goes beyond this by considering cases in which a non-linear mapping is applied to the output of the inner expectation, (e.g. the logarithm in the experimental design example), prohibiting such reformulation. We demonstrate that the construction of consistent NMC algorithms is possible, establish convergence rates, and provide empirical evidence that these rates are observed in practice. Our results show that whenever an outer estimator depends non-linearly on an inner estimator, then the number of samples used in both the inner and outer estimators must, in general, be driven to infinity for convergence. We extend our results to cases of repeated nesting and show that the optimal NMC convergence rate is $O(1/T^{\frac{2}{D+2}})$ where $T$ is the total number of samples used in the estimator and $D$ is the nesting depth (with $D=0$ being conventional MC), whereas naïve approaches only achieve a rate of $O(1/T^{\frac{1}{D+1}})$. We further lay out methods for reformulating certain classes of nested expectation problems into a single expectation, allowing usage of conventional MC estimation schemes with superior convergence rates than naïve NMC. Finally, we use our results to make application-specific advancements in Bayesian experimental design and variational auto-encoders.

Though the convergence of NMC has previously received little attention within the machine learning literature, a number of special cases having been investigated in other fields, sometimes under the name of nested simulation (Longstaff and Schwartz, 2001; Belomestny et al., 2010; Gordy and Juneja, 2010; Broadie et al., 2011). While most of this literature focuses on particular application-specific non-linear mappings, a convergence bound for a wider range of problems was shown by Hong and Juneja (2009) and recently revisited in the context of rare-event problems by Fort et al. (2017). The latter paper further considers the case where samples in the outer estimator originate from a Markov chain. Compared to this previous work, ours is the first to consider multiple levels of nesting, applies to a wider range of non-linear mappings, and provides more precise convergence rates. By introducing new results, outlining special cases, providing empirical assessment, and examining specific applications, we provide a unified investigation and practical guide nesting MC estimators in a machine learning context. We begin to realize the potential significance of this by using our theoretical results to make advancements in a number of specific application areas.

Another body of literature related to our work is in the study of the convergence of Markov chains with approximate transition kernels (Rudolf and Schweizer, 2015; Alquier et al., 2016; Medina-Aguayo et al., 2016). The analysis in this work is distinct, but complementary, to our own, focusing on the impact of a known bias on an MCMC chain, whereas our focus is more on the quantifying this bias. Also related is the study of techniques for variance reduction, such as multilevel MC (Heinrich, 2001; Giles, 2008), and bias reduction, such as the multi-step Richardson-Romberg method (Pages, 2007; Lemaire et al., 2017) and Russian roulette sampling (Lyne et al., 2015), many of which are applicable in a NMC context and can improve performance.

Problem Formulation

The key idea of MC is that the expectation of an arbitrary function $\lambda\colon\mathcal{Y}\rightarrow\mathcal{F}\subseteq\real$ under a probability distribution $p(y)$ for its input $y\in\mathcal{Y}$ can be approximated using:

In this paper, we consider the case that $\lambda$ is itself intractable, defined only in terms of a functional mapping of an expectation. Specifically, $\lambda(y)=f(y,\gamma(y))$ where we can evaluate $f\colon\mathcal{Y}\times\Phi\rightarrow\mathcal{F}$ exactly for a given $y$ and $\gamma(y)$, but $\gamma(y)$ is the output of the following intractable expectation of another variable $z\in\mathcal{Z}$:

depending on the problem, with $\phi\colon\mathcal{Y}\times\mathcal{Z}\rightarrow\Phi$. All our results apply to both cases, but we will focus on (3a) for clarity. Estimating $I$ involves computing an integral over $z$ for each value of $y$ in the outer integral. We refer to the approach of tackling both integrations using MC as nested Monte Carlo (NMC):

where each $z_{n,m}\sim p(z|y_{n})$ are independently sampled. In Section 3 we will build on this further by considering cases with multiple levels of nesting, where calculating $\phi(y,z)$ involves computation of an intractable (nested) expectation.

Convergence of Nested Monte Carlo

We now show that approximating $I\approx I_{N,M}$ is in principle possible, at least when $f$ is well-behaved. In particular, we establish a convergence rate of the mean squared error of $I_{N,M}$ and prove a form of almost sure convergence to $I$. We

further generalize our convergence rate to apply to the case of multiple levels of estimator nesting.

More formally, convergence bounds for NMC have previously been shown by Hong and Juneja (2009). Under the assumptions that each $\left(\hat{\gamma}_{M}\right)_{n}$ is Gaussian distributed (which is often reasonable due to the central limit theorem) and that $f$ is thrice differentiable other than at some finite number of points, they show that it is possible to achieve a converge rate of $O(1/N+1/M^{2})$. We now show that these assumptions can be relaxed to only requiring $f$ to be Lipschitz continuous, at the expense of weakening the bound.

If $f$ is Lipschitz continuous and $f(y_{n},\gamma(y_{n})),\phi(y_{n},z_{n,m})\in L^{2}$, the mean squared error of $I_{N,M}$ converges to at rate $O\left(1/N+1/M\right)$.

The theorem follows as a special case of Theorem 3. For exposition, a more accessible proof for this particular result is also provided in Appendix A in the supplement. ∎

Inspection of the convergence rate above shows that, given a total number of samples $T=MN$, our bound is tightest when $N\propto M$, with a corresponding rate $O(1/\sqrt{T})$ (see Appendix G). When the additional assumptions of Hong and Juneja (2009) apply, this rate can be lowered to $O(1/T^{2/3})$ by setting $N\propto M^{2}$. We will later show that this faster convergence rate can be achieved whenever $f$ is continuously differentiable, see also (Fort et al., 2017).

These convergence rates suggest that, for most $f$, it is necessary to increase not only the total number of samples, $T$, but also the number of samples used for each evaluation of the inner estimator, $M$, to achieve convergence. Further, as we show in Appendix B, the estimates produced by NMC are, in general, biased. This is perhaps easiest to see by noting that as $N\to\infty$, the variance of the estimator must tend to zero by the law of large numbers, but our bounds remain non-zero for any finite $M$, implying a bias.

2 Repeated Nesting and Exact Bounds

We next consider the case of multiple levels of nesting. As previously explained, this case is particularly important for analyzing probabilistic programming languages. To formalize what we mean by arbitrary nesting, we first assume some fixed integral depth $D>0$, and real-valued functions $f_{0},\cdots,f_{D}$. We then define

for $0\leq k\leq D-1$, where each $y^{(k)}_{n}\sim p\left(y^{(k)}|y^{(0:k-1)}\right)$ is drawn independently. Note that there are multiple values of $y^{(k)}_{n}$ for each possible $y^{(0:k-1)}$ and that $I_{k}\left(y^{(0:k-1)}\right)$ is still a random variable given $y^{(0:k-1)}$.

We are now ready to provide our general result for the convergence bounds that applies to cases of repeated nesting, provides constant factors (rather than just using big $O$ notation), and shows how the bound can be improved if the additional assumption of continuous differentiability holds.

If $f_{0},\cdots,f_{D}$ are all Lipschitz continuous in their second input with Lipschitz constants

where $O(\epsilon)$ represents asymptotically dominated terms.

If $f_{0},\cdots,f_{D}$ are also continuously differentiable with second derivative bounds

then this mean square error bound can be tightened to

For a single nesting, we can further characterize $O(\epsilon)$ giving

for when the continuous differentiability assumption does not hold and holds respectively.

These results give a convergence rate of $O(\sum_{k=0}^{D}1/N_{k})$ when only Lipschitz continuity holds and $O(1/N_{0}+(\sum_{k=1}^{D}1/N_{k})^{2})$ when all the $f_{k}$ are also continuously differentiable. As estimation requires drawing $O(T)$ samples where $T=\prod_{k=0}^{D}N_{k}$, the convergence rate will rapidly diminish with repeated nesting. More precisely, as shown in Appendix G, the optimal convergence rates are $O(1/T^{\frac{1}{D+1}})$ and $O(1/T^{\frac{2}{D+2}})$ respectively for the two cases, both of which imply that the rate diminishes exponentially with $D$.

Special Cases

We now outline some special cases where it is possible to achieve a convergence rate of $O(1/N)$ in the mean square error (MSE) as per conventional MC estimation. Establishing these cases is important because it identifies for which problems we can use conventional results, when we can achieve an improved convergence rate, and what precautions we must take to ensure this. We will focus on single nesting instances, but note that all results still apply to repeated nesting scenarios because they can be used to “collapse” layers and thereby reduce the depth of the nesting.

Our first special case is that $f$ is linear in its second argument, i.e. $f(y,\alpha v+\beta w)=\alpha f(y,v)+\beta f(y,w)$. Here the problem can be rearranged to a single expectation, a well-known result which forms the basis for pseudo-marginal, nested sequential MC (Naesseth et al., 2015), and certain ABC methods (Csilléry et al., 2010). Namely we have

where $(y_{n},z_{n})\sim p(y)p(z|y)$ if $\gamma(y)$ is of the form of (3a) and $y_{n}\sim p(y)$ and $z_{n}\sim p(z)$ are independently drawn if $\gamma(y)$ is of the form of (3b).

2 Finite Possible Realizations of y𝑦y

Our second case is if $y$ must take one of finitely many values $y_{1},\cdots,y_{C}$, then it is possible to use another approach to ensure the same convergence rate as standard MC. The key observation is to note that in this case we can convert the nested problem (2) into $C$ separate non-nested problems

with $y_{n}\overset{i.i.d.}{\sim}p(y)$ and $z_{n,c}\sim p(z|y_{c})$ (or $z_{n,c}\sim p(z)$ if using the formulation in (3b)). Note the critical point that each $z_{n,c}$ is independent of $y_{n}$ as each $y_{c}$ is a constant. We can now show the following result which, though intuitively straightforward, requires care to formally prove.

If $f$ is Lipschitz continuous, then the mean squared error of $I_{N}=\sum_{c=1}^{C}(\hat{P}_{N})_{c}\,(\hat{f}_{N})_{c}$ as an estimator for $I$ as per (10) converges at rate $O(1/N)$.

3 Products of Expectations

We next consider the scenario, which occurs for many latent variables models and probabilistic programming problems, where $\gamma(y)$ is equal to the product of multiple expectations, rather than just a single expectation as per (3a). That is,

which is a single expectation on an extended space and shows that (14) fits the NMC formulation. Furthermore, we can now show that if $f$ is linear, the MSE of the NMC estimator (14) converges at the standard MC rate $O(1/N)$, provided that $M$ remains fixed.

As this result holds in the case $L=1$, an important consequence is that whenever $f$ is linear, the same convergence rate is achieved regardless of whether we reformulate the problem to a single expectation or not, provided that the number of samples used by the inner estimator is fixed.

4 Polynomial f𝑓f

Perhaps surprisingly, whenever $f$ is of the form

where $z$ and $z^{\prime}$ are i.i.d. Therefore, assuming appropriate integrability requirements, we can construct the following non-nested MC estimator:

Empirical Verification

The convergence rates proven in Section 3 are only upper bounds on the worst-case performance. We will now examine whether these convergence rates are tight in practice, investigate what happens when our guidelines are not followed, and outline some applications of our results.

We start with the following analytically calculable problem

for which $I=\frac{1}{2}\log\left(\frac{2}{5\pi}\right)-\frac{2}{15}$. Figure 2(a) shows the corresponding empirical convergence obtained by applying (4) to (18) directly. It shows that, for this problem, the theoretical convergence rates from Theorem 3 are indeed realized. The figure also demonstrates the danger of not increasing $M$ with $N$, showing that the NMC estimator converges to an incorrect solution when $M$ is held constant. Figure 2(b) shows the effect of varying $N$ and $M$ for various fixed sample budgets $T$ and demonstrates that the asymptotically optimal strategy can be suboptimal for finite budgets.

2 Planning Cancer Treatment

3 Repeated Nesting

We next consider some simple models with multiple levels of nesting, starting with

which has analytic solution $I=-3/32$. The convergence plot shown in Figure 4 demonstrates that the theoretically expected convergence behaviors are observed for different methods of setting $N_{0},N_{1}$, and $N_{2}$.

As a byproduct, BOPP also produced Gaussian process approximations to the log MSE variations, as shown in Figure 5. We see that the two problems lead to distinct performance variations. Based on the (unshown) uncertainty estimates of these Gaussian processes, we believe these approximations are a close representation of the truth.

Applications

In this section, we show how our results can be used to derive an improved estimator for the problem of Bayesian experimental design (BED) in the case where the experiment outputs are discrete. A summary of our approach is provided here, with full details provided in Appendix I.

Bayesian experimental design provides a framework for designing experiments in a manner that is optimal from an information-theoretic viewpoint (Chaloner and Verdinelli, 1995; Sebastiani and Wynn, 2000). Given a prior $p(\theta)$ on parameters $\theta$ and a corresponding likelihood $p(y|\theta,d)$ for experiment outcomes $y$ given a design $d$, the Bayesian optimal design $d^{*}$ is given by maximizing the mutual information between $\theta$ and $y$ defined as follows

Estimating $d^{*}$ is challenging as $p(\theta|y,d)$ is rarely known in closed-form. However, appropriate algebraic manipulation shows that (21) is consistently estimated by

where $\theta_{n,m}\sim p(\theta)$ for each $(m,n)\in\{0,\ldots,M\}\times\{1,\ldots,N\}$, and $y_{n}\sim p(y|\theta=\theta_{n,0},d)$ for each $n\in\{1,\ldots,N\}$. This naïve NMC estimator has been implicitly used by (Myung et al., 2013) amongst others and gives a convergence rate of $O(1/N+1/M^{2})$ as per Theorem 3.

When $y$ can only take on finitely many realizations $y_{1},\dots,y_{c}$, we use the ideas introduced in Section 4.2 to derive the following improved estimator

where $\theta_{n}\sim p(\theta),\forall n\in\{1,\dots,N\}$. As $C$ is fixed, (23) converges at the standard MC error rate of $O(1/N)$. This constitutes a substantially faster convergence as (22) requires a total of $MN$ samples compared to $N$ for (23).

We finish by showing that the theoretical advantages of this reformulation also lead to empirical gains. For this we consider a model used in psychology experiments introduced by (Vincent, 2016), details of which are given in Appendix I. Figure 6 demonstrates that the theoretical convergence rates are observed while results given in Appendix I show that this leads to significant practical gains in estimating $\bar{U}(d)$.

2 Variational Autoencoders

To give another example of the applicability of our results, we now use Theorem 3 to directly derive a new result for the importance weighted autoencoder (IWAE) (Burda et al., 2015). Both the IWAE and the standard variational autoencoder (VAE) (Kingma and Welling, 2013) use lower bounds on the model evidence as objectives for train deep generative models and employ estimators of the form

3 Nesting Probabilistic Programs

Probabilistic programming systems (PPSs) (Goodman et al., 2008; Wood et al., 2014) provide a strong motivation for the study of NMC methods because many PPSs allow for arbitrary nesting of models (or queries, as they are known in the PPS literature), such that it is easy to define and run nested inference problems, including cases with multiple layers of nesting (Stuhlmüller and Goodman, 2012, 2014). Though this ability to nest queries has started to be exploited in application-specific work (Ouyang et al., 2016; Le et al., 2016), the resulting nested inference problems fall outside the scope of conventional convergence proofs and so the statistical validity of the underlying inference engines has previously been an open question in the field.

As we show in Rainforth (2017, 2018), the results presented here can be brought to bear on assessing the relative correctness of the different ways PPSs allow model nesting. In particular, the correctness of sampling from the conditional distribution of one query within another follows from Theorem 3, but only if the computation for each call to the inner query increases the more times that query is called. This requirement is not satisfied by current systems. Meanwhile, Theorem 5 can be used to the assert that observing the output of one query inside another leads to convergence at the standard MC rate, provided that the computation of the inner query instead remains fixed.

Conclusions

We have introduced a formal framework for NMC estimation and shown that it can be used to yield a consistent estimator for problems that cannot be tackled with conventional MC alone. We have derived convergence rates and considered what minimal continuity assumptions are required for convergence. However, we have also highlighted a number of potential pitfalls for naïve application of NMC and provided guidelines for avoiding these, e.g. highlighting the importance of increasing the number of samples in both the inner and the outer estimators to ensure convergence. We have further introduced techniques for converting certain classes of NMC problems to conventional MC ones, providing improved convergence rates. Our work has implications throughout machine learning and we hope it will provide the foundations for exploring this plethora of applications.

Appendix A Proof of Theorem 1 - Simplified Convergence Rate

Though the Theorem follows directly from Theorem 3, we also provide the following proof for this simplified case to provide a more accessible intuition behind the result. Note that the approach taken is distinct from the proof of Theorem 3.

Using Minkowski’s inequality, we can bound the mean squared error of $I_{N,M}$ by

We see immediately that $U=O\left(1/\sqrt{N}\right)$, since $\frac{1}{N}\sum_{n=1}^{N}f(y_{n},\gamma(y_{n}))$ is a MC estimator for $I$, noting our assumption that $f(y_{n},\gamma(y_{n}))\in L^{2}$. For the second term,

where $K$ is a fixed constant, again by Minkowski and using the assumption that $f$ is Lipschitz. We can rewrite

by the tower property of conditional expectation, and note that

since each $z_{n,m}$ is i.i.d. and conditionally independent given $y_{n}$. As such

Substituting these bounds for $U$ and $V$ in (25) gives

Appendix B The Inevitable Bias of Nested Estimation

In this section we demonstrate formally that NMC schemes must produce biased estimates of $I(f)$ for certain functions $f$. In fact, our result applies more generally: we show that this holds for any MC scheme that makes use of imperfect estimates $\hat{\zeta}_{n}$ of $\gamma(y_{n})$, either via a NMC procedure (e.g. $\hat{\zeta}_{n}=(\hat{\gamma}_{M})_{n}$), or when these inner estimates are generated by some other methods such as a variational approximation (Blei et al., 2016) or Bayesian quadrature (O’Hagan, 1991).

It also follows from Jensen’s inequality that any strictly convex or concave $f$ entails a biased estimator when $\hat{\zeta}_{n}$ is unbiased but has non-zero variance given $y_{n}$, e.g. when $\hat{\zeta}_{n}$ is a MC estimate. More formally we have

Similarly for any $f$ that is strictly concave in its second argument,

We prove our claim for the case that $f$ is strictly convex; our proof for the other concave case is symmetrical. We have

where the $\geq$ is a result of Jensen’s inequality on the inner expectation. Since $f$ is strictly convex and therefore non-linear, equality holds if and only if $\hat{\zeta}_{1}$ is almost surely constant given $y_{1}$. This is violated whenever $y_{1}\in\mathcal{A}$, which by assumption has a non-zero probability of occurring. Consequently, the inequality must be strict, giving the desired result. ∎

In addition to some special cases discussed in the Section 4, it may still be possible to develop unbiased estimation schemes for certain non-linear $f$ using Russian roulette sampling (Lyne et al., 2015) or other debiasing techniques. However, these induce their own complications: for some problems the resultant estimates have infinite variance (Lyne et al., 2015) and as shown by (Jacob et al., 2015), there are no general purpose “$f$-factories” that produce both non-negative and unbiased estimates for non-constant, positive output functions $f:\real\rightarrow{}^{+}$, given unbiased estimates for the inputs.

Appendix C Proof of Theorem 2 - “Almost almost sure” convergence

For all $N,M$, we have by the triangle inequality that

A second application of the triangle inequality then allows us to write

For any fixed $\delta>0$ then by repeatedly applying Egorov’s theorem to each $M\geq L$, we can find a sequence of events

for all $M\geq L$, $M^{\prime}\geq\tau^{1}_{\delta}(M)$, and $\omega\not\in B_{M}$, remembering that $\omega$ is a point in our sample space. We further have that (27) holds for all $M\geq M_{0}$, $M^{\prime}\geq\tau^{1}_{\delta}(M)$, and $\omega\not\in B_{\delta}:=\bigcup_{M\geq L}B_{M}$. Consequently, for all such $M$, $M^{\prime}$ and $\omega$,

To complete the proof, we must remove the dependence of $U_{N}$ on $N$ as well. This is straightforward once we observe that $U_{N}\overset{a.s.}{\to}0$ as $N\to\infty$ by the strong law of large numbers. So, by Egorov’s theorem again, there exists an event $C_{\delta}$ such that

We can now define $\tau_{\delta}(M)=\max(\tau^{1}_{\delta}(M),\tau^{2}_{\delta}(M))$, and $A_{\delta}=B_{\delta}\cup C_{\delta}$. By inequalities in (29) and (30),

Also, by the inequalities in (28) and (31),

Appendix D Proof of Theorem 3 - Convergence for Repeated Nesting

As this is a long and involved proof, we start by defining a number of useful terms that will be used throughout. Unless otherwise stated, these definitions hold for all $k\in\left\{0,\dots,D\right\}$. Note that most of these terms implicitly depend on the number of samples $N_{0},N_{1},\dots,N_{D}$. However, $s_{k}$, $\zeta_{d,k}$, and $\varsigma_{k}$ do not and are thus constants for a particular problem.

Given these definitions, we start by breaking the error down into a variance and bias term. Using the standard bias-variance decomposition we have

It is immediately clear from its definition in (36) that the bias term $\left(\beta_{k}\left(y^{(0:k-1)}\right)\right)^{2}$ is independent of $N_{0}$. On the other hand, we will show later that the dominant components of the variance term for large $N_{0:D}$ depend only on $N_{0}$. We can thus think of increasing $N_{0}$ as reducing the variance of the estimator and increasing $N_{1:D}$ as reducing the bias.

with the equality following because the $y_{n}^{(0:k)}$ being drawn i.i.d. and the expectation of each $f_{k}\left(y^{(0:k)},I_{k+1}\left(y^{(0:k)}\right)\right)$ equaling $\bar{f}_{k}\left(y^{(0:k-1)}\right)$ means that all the cross terms are zero. By the definition of $\sigma_{k}^{2}$ we now have

By using Minkowski’s inequality and the definition of $A_{k}$ it also follows that

Using a bias-variance decomposition on the second term above and noting that $s_{k}^{2}\left(y^{(0:k-1)}\right)$ and $\bar{f}_{k}\left(y^{(0:k-1)}\right)-\beta_{k}\left(y^{(0:k-1)}\right)$ are respectively the variance and expectation of $f_{k}\left(y^{(0:k)},\gamma_{k+1}\left(y^{(0:k)}\right)\right)$, we can rearrange the right-hand size of (47) to give

Here $s_{k}^{2}$ is independent of the number of samples used at any level of the estimate, while $A_{k}$ and $\beta_{k}^{2}$ are independent of $N_{d}\;\forall d\leq k$. Now by Jensen’s inequality, we have that

noting that the only difference in the definition of $\left(\beta_{k}\left(y^{(0:k-1)}\right)\right)^{2}$ and $A_{k}\left(y^{(0:k-1)}\right)$ is whether the squaring occurs inside or outside the expectation. Therefore, presuming that $A_{k}$ does not increase with $N_{d}\;\forall d>k$, neither will $\sigma_{k}^{2}\left(y^{(0:k-1)}\right)$, and so the variance term will converge to zero with rate $O(1/N_{k})$. Further, if ${A_{k}}\rightarrow 0$ as $N_{k+1},\dots,N_{D}\rightarrow\infty$, then for a large number of inner samples $\sigma_{k}^{2}\rightarrow s_{k}^{2}$ and thus we will have $v_{k}^{2}\left(y^{(0:k-1)}\right)\leq\frac{s_{k}^{2}}{N_{k}}+O\left(\epsilon\right)$ where $O\left(\epsilon\right)$ represents higher order terms that are dominated in the limit $N_{k},\dots,N_{D}\rightarrow\infty$. Provided this holds, we will also, therefore, have that

We now show that Lipschitz continuity is sufficient for ${A_{k}}\rightarrow 0$ and derive a concrete bound on the variance by bounding ${A_{k}}$. By definition of Lipschitz continuity, we have that

For the case where we only assume Lipschitz continuity then we will simply use the bound on the bias given by (49) leading to

which fully defines a bound on conditional the variance of one layer given the mean squared error of the layer below. In particular as $\omega_{D}\left(y^{(0:D-1)}\right)=0$ we now have

which is the standard error for Monte Carlo convergence. We further have

This leads to the following result for the single nesting case

$\approx\frac{\varsigma^{2}_{0}}{N_{0}}+\frac{K_{0}^{2}\varsigma_{1}^{2}}{N_{1}}=O\left(\frac{1}{N_{0}}+\frac{1}{N_{1}}\right)$ where the approximation becomes exact as $N_{0},N_{1}\rightarrow\infty$. Note that there is no $O\left(\epsilon\right)$ term as this bound is exact in the finite sample case.

Things quickly get messy for double nesting and beyond so we will ignore non-dominant terms in the limit $N_{0},\dots,N_{D}\rightarrow\infty$ and resort to using $O(\epsilon)$ for these instead. We first note that removing dominated terms from (53) gives

as $s_{k}^{2}$ does not decrease with increasing $N_{k+1:D}$ whereas the other terms do. We therefore also have

and therefore by recursively substituting (57) into itself we have

By definition we have that $\zeta_{0,0}^{2}=s_{0}^{2}=\varsigma_{0}^{2}$ and $\zeta_{d,0}^{2}=\varsigma_{d}^{2}$ and as (59) holds in the case $k=0$, the mean squared error of the overall estimator is as follows

and we have the desired result for the Lipschitz case.

We now revisit the bound for the bias in the continuously differentiable case to show that a tighter overall bound can be found. We first remember that

Taylor’s theorem implies that for any continuously differentiable $f_{k}$ we can write

where $\lim_{x\rightarrow\gamma_{k+1}\left(y^{(0:k)}\right)}h_{3}(x)=0$. Consequently, the last term is $O\left(\left(I_{k+1}\left(y^{(0:k)}\right)-\gamma_{k+1}\left(y^{(0:k)}\right)\right)^{3}\right)$ and so will diminish in magnitude faster than the first two terms provided that the derivatives are bounded, which is guaranteed by our assumptions. We will thus use $O(\epsilon)$ for this term and note that it is dominated in the limit.

By using the tower property we further have that

Remembering (50) we can recursively define the error bound in the same manner as the Lipschitz case. We can immediately see that, as $\beta_{D}=0$ without any nesting, we recover the bound from the Lipschitz case for the inner most estimator as expected. As the innermost estimator is unbiased we also have $\lambda_{D-1}\left(y^{(0:D-2)}\right)=0$ and so

Going back to our original bound on $\sigma_{D-1}^{2}\left(y^{(0:D-2)}\right)$ given in (48) and substituting for $\beta_{D-1}\left(y^{(0:D-2)}\right)$ we now have

There does not appear to be tighter bound for $A_{D-1}\left(y^{(0:D-2)}\right)$ than in the Lipschitz continuous case and so using the same bound of $A_{D-1}\left(y^{(0:D-2)}\right)\leq K_{D-1}^{2}\zeta^{2}_{D,D-1}\left(y^{(0:D-2)}\right)/N_{D-1}$ we have

Therefore for the single nesting case, we now have

$\approx\frac{\varsigma^{2}_{0}}{N_{0}}+\frac{C_{0}^{2}\varsigma_{1}^{4}}{4N_{1}^{2}}=O\left(\frac{1}{N_{0}}+\frac{1}{N_{1}^{2}}\right)$ where again the approximation becomes tight when $N_{0},N_{1}\rightarrow\infty$. Here we have used the fact that the only $O(\epsilon)$ term comes from the Taylor expansion and is equal to $O\left(\frac{1}{N_{1}^{3}}\right)$ because we have $\delta_{1,D-1}=0$ and therefore

Returning to calculating the bound for the repeated nesting case then by substituting (62) into (50) we have more generally

We also have that except at $k=D-1$ and $k=D$ (for which both $\lambda_{k}$ and $\beta_{k+1}$ are zero), then

for sufficiently large $N_{k+1},\dots,N_{D}$. This means that when we substitute (67) into (66), the second term in (67) becomes dominated giving

Now as $\beta_{k+1}^{2}\left(y^{(0:k)}\right)=E_{k+1}\left(y^{(0:k)}\right)-\frac{s_{k+1}^{2}\left(y^{(0:k)}\right)}{N_{k+1}}$ we have

Appendix E Proof of Theorem 4 - Convergence Rate for Finite Realisations of y𝑦y

Denote $P_{c}=P(y=y_{c})$ and $f_{c}=f(y_{c},\gamma(y_{c}))$ noting that as the $y_{c}$ are fixed values, so are $P_{c}$ and $f_{c}$. Then, Minkowski’s inequality allows us to bound the mean squared error as

Moreover, again by Minkowski, we have $W_{c}\leq U_{c}+V_{c}$ where

Factoring out $(\hat{P}_{N})_{c}$ in $U_{c}$ and noting that each $y_{n}$ and $z_{n,c}$ are sampled independently gives

Using Minkowski’s inequality, we may write the first right-hand term as

For the second term, note that by Lipschitz continuity, we have for some constant $K>0$

since $\frac{1}{N}\sum_{n=1}^{N}\phi(y_{c},z_{n,c})$ is a Monte Carlo estimator for $\gamma(y_{c})$. Altogether then, we have that

We can also factor out $f_{c}$ in $V_{c}$ to obtain

Appendix F Proof for Theorem 5 - Products of Expectations

Appendix G Optimizing the Convergence Rates

We have shown that the mean squared error converges at a rate

depending on the smoothness assumptions that can be made about $f$. Here we show that given a sample budget for the inner most estimator $T=\prod_{k=0}^{D}N_{k}$, then these bounds are optimized by setting $N_{0}\propto N_{1}\propto\dots\propto N_{D}$ and $N_{0}\propto N_{1}^{2}\propto\dots\propto N_{D}^{2}$ respectively for the two cases and that this gives bounds of $O\left(1/T^{\frac{1}{D+1}}\right)$ and $O\left(1/T^{\frac{2}{D+2}}\right)$respectively. For the single nested case, this gives bounds of $O(1/\sqrt{T})$ and $O(1/T^{2/3})$ respectively.

We start by explaining why $T$ is an appropriate measure of the overall computational cost. First note that for each sample of $y^{(0:k)}$, the NMC estimator requires $N_{k}$ samples of $y^{(k+1)}$. Thus there are $N_{0}$ samples of the outermost level, $N_{0}\times N_{1}$ of the next level, and $T=\prod_{k=0}^{D}N_{k}$ samples of the innermost level, regardless of the setup. In other words, each individual estimate of the innermost level uses $N_{D}$ samples and we generate $\prod_{k=0}^{D-1}N_{k}=T/N_{D}$ of these estimates because we need to generate one estimate for each sample of the layer above. Thus what we can vary for a fixed $T$ is whether we use more estimates each using fewer samples, or fewer estimates each using more samples.

Now the total cost of generating $I_{0}$ scales with sum the costs of each individual layer, namely

To derive the optimal rates, we first consider the single nested case where $D=1$, $N_{0}=N$, and $N_{1}=M$. Consider setting $N=\tau(M)$ then $T=\tau(M)\cdot M$ and our bounds become $O(R)$, where

In this first case supposing $\tau(M)=O(M)$ easily gives

as $M\to\infty$. In contrast, consider the case that $\tau(M)\gg M$ as $M\to\infty$. We then have $\frac{1}{\sqrt{M}}\gg\frac{1}{\sqrt{\tau(M)}}$ as $M\to\infty$, so that

Comparing with (69), we observe that, for the same total budget of samples $T$, this choice of $\tau$ provides a strictly weaker convergence guarantee than in the previous case. When $M\gg\tau(M)$ also then we have $\frac{1}{\sqrt{\tau(M)}}\gg\frac{1}{\sqrt{M}}$ as $M\to\infty$ and so

which is again a weaker bound. We thus see that the $O(1/N+1/M)$ bound is optimized when $N\propto M$, giving a convergence rate of $O(1/\sqrt{T})$.

In the second case suppose that $\tau(M)=O(M^{2})$ as $M\to\infty$. This now gives

as $M\to\infty$. Now considering the cases $\tau(M)\gg M^{2}$ leads to $\frac{1}{M^{4/3}}\gg\frac{1}{\tau(M)^{2/3}}$ and thus

Similarly, if $\tau(M)\ll M^{2}$ then $\frac{1}{\tau(M)^{1/3}}\gg\frac{1}{M^{2/3}}$ and thus

Both of these cases lead to weaker bounds and so we see that the $O(1/N+1/M^{2})$ bound is tightest when $N\propto M^{2}$, giving a convergence rate of $O(1/T^{2/3})$.

We now consider the repeated nesting case without continuously differentiability such that our bound is $O\left(\sum_{k=0}^{D}\frac{1}{N_{k}}\right)$. Here we can immediately see that $N_{0}\propto N_{1}\propto\dots\propto N_{D}$ leads to $N_{k}\propto T^{\frac{1}{D+1}}$ and thus $O\left(1/T^{\frac{1}{D+1}}\right)$ convergence. If we were to set any $N_{k}\ll T^{\frac{1}{D+1}}$ then this term would dominate the sum and lead to a worse converge. Thus the result from the single nested case trivially extends to the multiple nested case, giving the required result.

Finally considering repeated nesting for the bound $O\left(\frac{1}{N_{0}}+\left(\sum_{k=1}^{D}\frac{1}{N_{k}}\right)^{2}\right)$ then we have from the previous result that $N_{1}\propto N_{2}\propto\dots\propto N_{D}$ is required for optimality as otherwise one of the terms in the summation dominates the other terms. If we now define $M=\prod_{k=1}^{D}N_{k}=T/N_{0}$ then we get a convergence rate of $O(1/N_{0}+1/M^{2})$ which is identical to the single nesting case for this tighter bound. We, therefore, have that the optimal configuration must be $N_{0}\propto N_{1}^{2}\propto\dots\propto N_{D}^{2}$ giving a bound of $O\left(1/T^{\frac{2}{D+2}}\right)$ as it gives $N_{0}\propto T^{\frac{2}{D+2}}$.

Appendix H Additional details pertaining to cancer simulator

In this section, we elucidate some more details about the cancer simulator described in the manuscript, provide more rigorous mathematical definitions for the relevant terms using the same nomenclature, and also include more results figures.

We define $I(T_{\text{treat}})$ to be the expected proportion of patients who receive treatment. A particular patient is represented by $y\in\mathcal{R}^{d}$. Specifically, $y$ consists of only a single real number ($d=1$) representing the size of the tumor upon discovery. Initial tumor size is drawn from a scaled Rayleigh distribution. The outcome of the simulator is then $\phi(y,z)\in\left\{0,1\right\}$, and is the binary outcome of whether that particular patient and sample of unobserved parameters yield an expected tumor size below the threshold, $T_{opp}$, after a fixed time duration, $t_{max}$. The simulator is a pair of coupled, parameterized differential equations for the action of an anti-tumor treatment such as chemotherapy, as described in Enderling and Chaplain (2014):

where $c(t,x)\in\mathcal{R}_{+}$ represents tumor size, with initial size $y_{n}$. Similarly, $K(t,x)\in\mathcal{R}_{+}$ represents the notion of a carrying capacity, with the initial carrying capacity, $K(0,z)$, set to a known constant $K_{0}$. The magnitude of the patient response to an anti-tumor treatment (such as chemotherapy) is represented by $\xi\in$, drawn from a beta distribution. $\{\lambda,\psi,\phi\}\in\mathcal{R}_{+}^{3}$ represent the parameters of the simulator. We also define $x_{n,m}=\{\lambda,\psi,\phi\,K_{0},\xi\}$ and $z_{n,m}=\{x_{n,m},T_{\text{opp}},t_{\text{max}}\}$, where all but £$\xi$ are set to constant values. Expanding this to condition all values on $y_{n}$ is trivial given domain knowledge. Alternatively, they could also be drawn at random, but not be conditioned on $y_{n}$. Such relations are omitted here for simplicity.

Taking the expectation of $\phi$ over $M$ different realizations of $z$ yields the estimate $(\hat{\gamma}_{M})_{n}$. This value is the probability that treatment will be successful for a particular patient, marginalizing over possible unobserved dynamics. This is the point at which clinician decides whether initiate the treatment plan. This decision is represented $f(y_{n},(\hat{\gamma}_{M})_{n})\in$ as:

where $T_{\text{treat}}$ is the minimum probability of success required for that patient to receive the treatment, and again, could be conditioned on $y$ also. Taking the expectation of $f$ over patients yields the expected frequency with which the treatment will be delivered, given a value of $T_{\text{treat}}$. The hospital wishes to estimate the value $T_{\text{treat}}$ that maximizes the number of patients treated, while only treating those patients with the highest probability of success, and (in expectation) staying within the budgetary constraint.

The model is completed by the definition of the following distributions and parameters.

H.2 Budget result

In the example outlined above, the treatment center is not actually attempting to evaluate the value of $I$, but to find the optimal value of $T_{\text{treat}}$ subject to a budgetary constraint. A simplistic way of evaluating the optimal value is to perform a dense search over different values of the parameter, each time evaluating the estimated expenditure, and select the best performing value.

Figure 7 shows the variation of predicted expenditure against the threshold probability, as well as the budget constraint. The intersection of these curves is the optimal setting of $T_{opp}$, here evaluated to be 12.5%. From the blue line, it is clear that the relationship between expenditure and treatment probability is non-linear, especially at the extrema of the distribution, and hence the use of NMC was necessarily for evaluating the optimal value.

Appendix I Bayesian Experimental Design

Bayesian experimental design provides a framework for designing experiments in a manner that is optimal from an information-theoretic viewpoint (Chaloner and Verdinelli, 1995; Sebastiani and Wynn, 2000). By minimizing the entropy in the posterior distribution of the parameters of interest, one can maximize the information gathered by the experiment.

Let the parameters of interest be denoted by $\theta\in\Theta$ for which we define a prior distribution $p(\theta)$. Let the probability of achieving outcome $y\in\mathcal{Y}$, given parameters $\theta$ and a design $d\in\mathcal{D}$, be defined by likelihood model $p(y|\theta,d)$. Under our model, the outcome of the experiment given a chosen $d$ is distributed according to

where we have used the fact that $p(\theta)=p(\theta|d)$ because $\theta$ is independent of the design. Our aim is to choose the optimal design $d$ under some criterion. We, therefore, define a utility function, $U(y,d)$, representing the utility of choosing a design $d$ and getting a response $y$. Typically our aim is to maximize information gathered from the experiment, and so we set $U(y,d)$ to be the gain in Shannon information between the prior and the posterior:

However, we are still uncertain about the outcome. Thus, we use the expectation of $U(y,d)$ with respect to $p(y|d)$ as our target:

noting that this corresponds to the mutual information between the parameters $\theta$ and the observations $y$. The Bayesian-optimal design is then given by

Finding $d^{*}$ is challenging because the posterior $p(\theta|y,d)$ is rarely known in closed form. To solve the problem, we proceed by rearranging (76) using Bayes’ rule (remembering that $p(\theta)=p(\theta|d)$):

The first of these terms can now be evaluated using standard MC approaches as the integrand is analytic. In contrast, the second term is not directly amenable to standard MC estimation as the marginal $p(y|d)$ represents an expectation and taking its logarithm represents a non-linear functional mapping.

To derive an estimator, we will now consider these terms separately. Starting with the first term,

where $\theta_{n}\sim p(\theta)$ and $y_{n}\sim p(y|\theta=\theta_{n},d)$. We note that evaluating (79) involves both sampling from $p(y|\theta,d)$ and directly evaluating it point-wise. The latter of these cannot be avoided, but in the scenario where we do not have direct access to a sampler for $p(y|\theta,d)$, we can use the standard importance sampling trick, sampling instead $y_{n}\sim q(y|\theta=\theta_{n},d)$ and weighting the samples in (79) by $w_{n}=\frac{p(y_{n}|\theta_{n},d)}{q(y_{n}|\theta_{n},d)}$.

where $\theta_{n,m}\sim p(\theta)$ and $y_{n}\sim p(y|d)$. Here we can sample the latter by first sampling an otherwise unused $\theta_{n,0}\sim p(\theta)$ and then sampling $y_{n}\sim p(y|\theta_{n,0},d)$. Again we can use importance sampling if we do not have direct access to a sampler for $p(y|\theta_{n,0},d)$.

Putting (79) and (80) together (and renaming $\theta_{n}$ from (79) as $\theta_{n,0}$ for notational consistency with (80)) we now have the following complete estimator given in the main paper and implicitly used by (Myung et al., 2013) amongst others

where $\theta_{n,m}\sim p(\theta)\;\forall m\in 0:M,\;n\in 1:N$ and $y_{n}\sim p(y|\theta=\theta_{n,0},d)\;\forall n\in 1:N$.

We now show that if $y$ can only take on one of $C$ possible values ($y_{1},\ldots,y_{C}$), we can achieve significant improvements in the convergence rate by using a similar to that introduced in Section 3.2 to convert to single MC estimator:

where $\theta_{n}\sim p(\theta)\quad\forall n\in 1,\dots,N$. As $C$ is a fixed constant, the MSE for first term clearly converges at the standard MC error rate of $O(1/N)$. Similarly each $\hat{P}_{N}(y_{c}|d)=\frac{1}{N}\sum_{n=1}^{N}p(y_{c}|\theta_{n},d)$ term also converges at a rate $O(1/N)$ to $p(y_{c}|d)$. Now noting that $\hat{P}_{N}(y_{c}|d)\leq 1$ and that $f(x)=x\log x$ is Lipschitz continuous in the range $(0,1]$, each $\hat{P}_{N}(y_{c}|d)\log\left(\hat{P}_{N}(y_{c}|d)\right)$ term must also converge at the MC error rate if $p(y_{c}|d)>0\;\forall c=1,\dots,C$. Finally if we assume that when $p(y_{c}|d)=0$ then $\hat{P}_{N}(y_{c}|d)=0$ almost surely for sufficiently large $N$, then the second term also converges at the MC error when $p(y_{c}|d)=0$. We now have a finite sum of terms which each convergence to $\bar{U}(d)$ with MC MSE rate $O(1/N)$, and so the overall estimator (82) must also converge at this rate. This compares to $O(1/T^{2/3})$ for (81) (assuming we take $N\propto M^{2}$), noting that generating $T$ samples for (81) has the same cost up to a constant factor as generating $N$ for (82). To the best of our knowledge, this is the first introduction of this superior estimator in the literature.

We finish by showing that the theoretical advantages of this reformulation also leads to empirical gains in the estimation of $\bar{U}(d)$. For this, we consider a model used in psychology experiments for delay discounting introduced by (Vincent, 2016; Vincent and Rainforth, 2017). Our experiment comprises of asking questions of the form “Would you prefer $\pounds A$ now, or $\pounds B$ in $D$ days?” and we wish to choose the question variables $d=\{A,B,D\}$ in the manner that will give the most incisive questions. The target participant is presumed to have parameters $\theta=\{k,\alpha\}$ and the following response model

where $y=1$ indicates choosing the delayed response and $\Phi$ represents the cumulative normal distribution. As more questions are asked, the distribution over the parameters $\theta$ is updated, such that the most optimal question to ask at a particular time depends on the previous questions and responses. For the sake of brevity, when comparing the performance of (81) and (82) we will neglect the problem of how best to optimize the design, and consider only the problem of evaluating $\bar{U}(d)$. We will further consider the case where $B=100$ and $D=50$ are fixed and we are only choosing the delayed value $A$. We presume the following distribution on the parameters

We first consider convergence in the estimate of $\bar{U}(d)$ for the case $A=70$ for our suggested method (82) and the naïve solution (81), the results of which are shown in Figure 2a in the main paper. Here we see that the convergence rates of the two methods are both as expected and that our suggested method offers significant empirical performance improvements.

We next consider setting a total sample budget $T=10^{4}$ and look at the variation in the estimated values of $\bar{U}(d)$ for different values of $A$ for the two methods as shown in Figure 8. This shows that the improvement in MSE leads to clearly visible improvements in the characterization of $\bar{U}(d)$ that will translate to improvements in seeking the optimum.

Acknowledgements

Tom Rainforth’s research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) ERC grant agreement no. 617071. However, the majority of this work was undertaken while he was in the Department of Engineering Science, University of Oxford, and was supported by a BP industrial grant. Robert Cornish is supported by an NVIDIA scholarship. Hongseok Yang is supported by an Institute for Information & communications Technology Promotion (IITP) grant funded by the Korea government (MSIP) (No.R0190-16-2011, Development of Vulnerability Discovery Technologies for IoT Software Security). Frank Wood is supported under DARPA PPAML through the U.S. AFRL under Cooperative Agreement FA8750-14-2-0006, Sub Award number 61160290-111668.

References