# Optimal randomized multilevel Monte Carlo for repeatedly nested expectations

Yasa Syed
Department of Statistics
Rutgers University  Guanyang Wang
Department of Statistics
Rutgers University

###### Abstract

The estimation of repeatedly nested expectations is a challenging task that
arises in many real-world systems. However, existing methods generally suffer
from high computational costs when the number of nestings becomes large. Fix
any non-negative integer $D$ for the total number of nestings. Standard Monte
Carlo methods typically cost at least ${\mathcal{O}}(\varepsilon^{-(2+D)})$
and sometimes ${\mathcal{O}}(\varepsilon^{-2(1+D)})$ to obtain an estimator up
to $\varepsilon$-error. More advanced methods, such as multilevel Monte Carlo,
currently only exist for $D=1$. In this paper, we propose a novel Monte Carlo
estimator called $\mathsf{READ}$, which stands for “Recursive Estimator for
Arbitrary Depth.” Our estimator has an optimal computational cost of
${\mathcal{O}}(\varepsilon^{-2})$ for every fixed $D$ under suitable
assumptions, and a nearly optimal computational cost of
${\mathcal{O}}(\varepsilon^{-2(1+\delta)})$ for any $0<\delta<\frac{1}{2}$
under much more general assumptions. Our estimator is also unbiased, which
makes it easy to parallelize. The key ingredients in our construction are an
observation of the problem’s recursive structure and the recursive use of the
randomized multilevel Monte Carlo method.

Keywords: nested expectation, optimal cost, randomized Multilevel Monte Carlo,
unbiased estimator

##  1 Introduction

Monte Carlo methods are a class of algorithms that use random sampling to
estimate quantities of interest, such as integrals or expected values. When
the estimand can be expressed as an expectation, for example
$\mathbf{E}_{\pi}[g(X)]$, these methods work by generating independent random
samples $X_{1},\ldots,X_{n}$ from $\pi$, and using the average
$\sum_{i=1}^{n}g(X_{i})/n$ as an estimator. Monte Carlo estimators are
unbiased and converge at a rate of $n^{-1/2}$, regardless of the dimension of
the samples. This dimension-independent convergence rate makes Monte Carlo
methods a powerful tool for approximating high-dimensional integrations, as
they do not suffer from the curse of dimensionality that plagues deterministic
numeric integration methods.

However, the above analysis implicitly assumes the integrand $g$ can be
pointwisely evaluated, which may not be possible in many situations of
practical interest. In this paper, we study the problem of estimating
repeatedly nested expectations (RNE), which means the integrand depends on a
sequence of other functions and conditional expectations. Specifically, fix
any positive integer $D$ for the total number of nestings, and
$\\{g_{d}\\}_{d=0}^{D}$ for a family of real-valued functions which can be
pointwisely evaluated. Let $(y^{(0)},\ldots,y^{(D)})$ be a finite-time
stochastic process with underlying joint distribution $\pi$, and let
$y^{(0:d)}$ denote the vector $(y^{(0)},\ldots,y^{(d)})$ for every $d\leq D$.
The RNE, first formally formulated in Rainforth et al. [2018], is defined as:

$\displaystyle\gamma_{0}=\mathbf{E}\left[g_{0}\left(y^{(0)},\gamma_{1}\left(y^{(0)}\right)\right)\right],$
(1)

where $\\{\gamma_{i}\\}_{i=1}^{D-1}$ is recursively defined as:

$\displaystyle\gamma_{d}(y^{(0:d-1)})=\mathbf{E}\left[g_{d}\left(y^{(0:d)},\gamma_{d+1}\left(y^{(0:d)}\right)\right)\bigg{|}~{}y^{(0:d-1)}\right],$
(2)

and

$\displaystyle\gamma_{D}(y^{(0:D-1)})=\mathbf{E}\left[g_{D}\left(y^{(0:D)}\right)\bigg{|}~{}y^{(0:D-1)}\right].$
(3)

Estimating RNEs is a fundamental problem that covers a variety of real-world
applications, where the quantity of interest depends on multiple stages or
decision points. For example, when
$g_{d}(y^{(0:d)},u)\coloneqq\max\\{y^{(d)},u\\}$ for $0\leq d\leq D-1$ and
$g_{D}(y^{(0:D)})=y^{(D)}$, the quantity $\gamma_{0}$ stands for the expected
utility of the optimal strategy in a $D$-horizon optimal stopping problem – a
central problem in financial modeling. Other applications include
probabilistic programs Rainforth [2018], Bayesian experimental design Goda et
al. [2022], portfolio risk management Gordy and Juneja [2010], physics and
chemistry Dauchet et al. [2018].

However, estimating RNEs is challenging. As shown in formulas (1) – (3), we
are interested in the expectation of $g_{0}$, which depends on the random
variable $y^{(0)}$ and $\gamma_{1}(y^{(0)})$ – a conditional expectation of
$g_{1}$ given $y^{(0)}$. Then $g_{1}$ further depends on a random variable
$y^{(1)}$ and $\gamma_{2}(y^{(0)},y^{(1)})$ which is a conditional expectation
of $g_{2}$ given $y^{(0)}$ and $y^{(1)}$. This procedure is recursively
defined until it reaches the deepest depth, $D$. Since $\gamma_{1}(y^{(0)})$
(and also $\gamma_{2},\gamma_{3},\ldots$) cannot be directly evaluated in most
practical cases, estimating RNEs cannot be handled by standard Monte Carlo
methods.

The most natural way to estimate RNEs is by nesting Monte Carlo (NMC)
estimators. In the $D=1$ case, this method works by first sampling independent
and identically distributed (i.i.d.) copies
$y^{(0)}_{1},\ldots,y^{(0)}_{N_{0}}$ according to the distribution of
$y^{(0)}$. For each fixed $y^{(0)}_{i}$, one further samples $N_{1}$ i.i.d.
$y^{(1)}_{1},\ldots,y^{(1)}_{N_{1}}$ according to $\pi(y^{(1)}\mid
y^{(0)}_{i})$, and uses the standard estimator
$\hat{\gamma}_{1}(y^{(0)}_{i})\coloneqq\sum_{i=j}^{N_{1}}g_{1}(y^{(0)}_{i},y^{(1)}_{j})/N_{1}$
to estimate $\gamma_{1}(y^{(0)}_{i})$. The final estimator uses the estimated
$\hat{\gamma}_{1}(y^{(0)}_{i})$ to replace the intractable
$\gamma_{1}(y^{(0)}_{i})$, i.e.,

$I_{N_{0},N_{1}}=\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}g_{0}(y^{(0)}_{i},\hat{\gamma_{1}}(y^{(0)}_{i})).$

This nested estimator can be easily extended to the general $D$ case, albeit
the notations become more complex. Roughly, one still samples $N_{0}$ i.i.d.
copies according to $\pi(y^{(0)})$, and for each fixed trajectory
$y^{(0:d-1)}$, the user generates $N_{d}$ i.i.d. samples from $\pi(y^{(d)}\mid
y^{(0:d-1)})$ all the way to depth $D$ and then form the nested estimator from
the deepest depth to the shallower depths. The construction details are
referred to Section 3.2 of Rainforth et al. [2018].

After suitably allocating the number of samples $(N_{d})_{d=0}^{D}$ for each
depth, the root-mean-square error (rMSE) of the NMC estimator converges to $0$
at a rate of $N^{-1/(2D+2)}$ or $N^{-1/(D+2)}$ Rainforth et al. [2018],
depending on the regularity conditions of the functions
$\\{g_{d}\\}_{d=0}^{D}$, where $N=\prod_{d=0}^{D}N_{d}$ is the total number of
samples used to form a nested estimator. This convergence rate diminishes
exponentially with $D$, meaning that NMC estimators do not have the same
dimension-free convergence rate as standard Monte Carlo estimators. As a
result, NMC methods require at least ${\mathcal{O}}(\varepsilon^{-(2+D)})$ and
sometimes ${\mathcal{O}}(\varepsilon^{-2(1+D)})$ samples to get an estimator
within $\varepsilon$-rMSE, while standard Monte Carlo estimators require only
$O(\varepsilon^{-2})$ samples. Although there are a few cases mentioned in
Rainforth et al. [2018] where the canonical ${\mathcal{O}}(N^{-1/2})$ rate can
be achieved, the problem of estimating RNEs with an optimal (or dimension-
free) convergence rate remains largely open.

In the special case $D=1$, more efficient methods have been proposed Giles
[2018], Giles and Goda [2019], Giles and Haji-Ali [2019] based on the
celebrated multilevel Monte Carlo (MLMC) methods Heinrich [2001], Giles
[2008]. These estimators achieve up to $\varepsilon$-rMSE with cost
$O(\varepsilon^{-2}\log(1/\varepsilon)^{2})$ or $O(\varepsilon^{-2})$ under
varying conditions, comparing favorably with the NMC estimator. However,
existing methods cannot be directly generalized to solve the general $D$ case.
Meanwhile, implementing these methods requires users to prespecify the
precision level $\varepsilon$ and conduct preliminary experiments/calculations
to carefully estimate/bound the parameters in the MLMC algorithm (see, e.g.,
Theorem 1 of Giles and Goda [2019]). Therefore, existing MLMC estimators seem
to be harder to implement and less amendable to our original problem, which
has a recursive structure.

In this work, we propose the $\mathsf{READ}$, a novel Monte Carlo estimator
for the RNE estimation with an arbitrary number of nestings $D$. Our
construction is interesting in the following three aspects. Firstly, under
suitable regularity conditions similar to those in Rainforth et al. [2018],
the rMSE of our estimator has an optimal convergence rate $N^{-1/2}$
regardless of $D$. Equivalently, our method costs in expectation
$O(\varepsilon^{-2})$ to get an estimator up to $\varepsilon$-rMSE. Under much
more general assumptions, our method still achieves a nearly-optimal cost of
${\mathcal{O}}(\varepsilon^{-2(1+\delta)})$ for any $0<\delta<\frac{1}{2}$ to
get an estimator up to $\varepsilon$-mean-absolute-error (MAE).

It is worth mentioning that most of our effort is devoted to designing
unbiased estimators of $\gamma_{0}$ in (1) with finite computational cost and
finite variance (or finite (2-$\delta$)-th moment under more general
assumptions). After developing such an unbiased estimator, we can simulate
independent copies of these estimators and average them. The $N^{-1/2}$
convergence rate and $\mathcal{O}(\varepsilon^{-2})$ cost are then immediate
corollaries of the bias-variance decomposition formula, see Corollary 2.3.

Therefore, another appealing property of $\mathsf{READ}$, in contrast to
existing methods, is that it admits no estimation bias. Unbiasedness implies
these estimators can be implemented in parallel processors without requiring
any communication between them. Designing unbiased estimators has recently
attracted much interest in statistics, operations research, and machine
learning communities for its potential for parallelization. Our methods add to
the rich body of works of Glynn and Rhee [2014], Rhee and Glynn [2015],
Blanchet and Glynn [2015], Jacob et al. [2020], Biswas et al. [2019], Wang et
al. [2021], Wang and Wang [2022], Kahale [2022].

Finally, our algorithm for constructing $\mathsf{READ}$ relies on the
randomized multilevel Monte Carlo (rMLMC) method McLeish [2011], Rhee and
Glynn [2015], Blanchet and Glynn [2015], but it is significantly different
from its previous applications. Most existing rMLMC applications Rhee and
Glynn [2015], Vihola [2018], Goda et al. [2022] also have a non-randomized
version (i.e., the original MLMC Giles [2008]) with similar or better
computational cost guarantees, leading some to believe that every problem
solved by rMLMC also has a natural non-randomized counterpart. However, our
work seems to suggest that this belief may not always hold true. The rMLMC
framework is well-suited to the recursive structure of RNEs, and can be used
as a subroutine in our method. In contrast, the non-randomized MLMC cannot be
easily applied to the general case of $D>1$. This suggests that the rMLMC
framework may be more widely applicable than previously thought.

The rest of this paper is organized as follows: in the remainder of this
section, we discuss related works, set up our notation, and introduce our
technical assumptions. In Section 2, we introduce our algorithm and show that
it attains the optimal and nearly-optimal computational cost under two
different assumptions, respectively. In Section 3, we demonstrate the
empirical performance of our method on several toy examples. We conclude this
paper with a short discussion in Section 4. Proof and experiment details are
deferred to the Appendix. An additional experiment is also included in
Appendix F.

###  1.1 Related work

Our algorithm design strategy mainly follows the randomized multilevel Monte
Carlo (rMLMC) framework McLeish [2011], Rhee and Glynn [2015], Blanchet and
Glynn [2015]. Our algorithm is inspired by the unbiased optimal stopping
estimator Zhou et al. [2022], which develops estimators for the optimal
stopping problem by recursively calling the rMLMC algorithm. We extend the
methodology in Zhou et al. [2022] both in scope and depth. Our method works
with a more general class of problems formulated by Rainforth et al. [2018],
which includes the optimal stopping problem as a special case, and provides
more precise results under practical assumptions.

Throughout this paper, we will assume the functions $\\{g_{d}\\}_{d=0}^{D-1}$
are all continuous and the process $\pi$ can be perfectly simulated. When
$D=1$ and $g_{0}$ is discontinuous, progress has been made by Broadie et al.
[2011] and Giles and Haji-Ali [2019, 2022]. When the underlying distribution
is itself challenging, users have to first use MCMC to approximately sample
from $\pi$. The case of $D=1$ and challenging $\pi$ is considered in Wang and
Wang [2022].

###  1.2 Notations

Now we introduce our notations. Many of our notations follow those used in the
original definition Rainforth et al. [2018]. Throughout this paper, we
preserve the letter $D$ for the total number of nestings. We denote by $\pi$
the underlying joint distribution of a finite-time, real-valued stochastic
process $(y^{(0)},\ldots,y^{(D)})$. For every $0\leq i\leq j\leq D$, we use
the $y^{(i:j)}$ to denote the vector $(y^{(i)},\ldots,y^{(j)})$. The
conditional distribution of $y^{(d:D)}$ given the value of $y^{(0:d-1)}$ is
denoted by $\pi_{d:D}(\cdot\mid y^{(0:d-1)})$. The marginal distribution of
$y^{(d)}$ conditioning on $y^{(0:d-1)}$ is denoted by $\pi_{d}(\cdot\mid
y^{(0:d-1)})$. We adopt the convention that $y^{(0:-1)}=\varnothing$, and
therefore $\pi_{0}$ stands for the (unconditioned) marginal distribution of
$y^{(0)}$. Let $\Pi$ be any probability distribution on some probability
space, and $Z$ be some random variable on the same space, then we use $\lVert
Z\rVert_{\Pi,m}$ to denote the $L^{m}$–norm of $Z$ under $\Pi$, i.e.,
$\big{(}\mathbf{E}_{\Pi}[\lvert Z\rvert^{m}]\big{)}^{1/m}$. The geometric
distribution with parameter $r$ is denoted by ${\mathsf{Geo}}(r)$. We also
define $p_{r}(n)\coloneqq\mathbf{P}[{\mathsf{Geo}}(r)=n]=r(1-r)^{n}$ for every
$n\in\\{0,1,2,\ldots,\\}$. For every $0\leq d<D$, the function $g_{d}$
introduced in (1) – (2) maps from $\mathbf{R}^{d+2}$ to $\mathbf{R}$. The
function $g_{D}$ in (3) maps from $\mathbf{R}^{D+1}$ to $\mathbf{R}$. For
$i.i.d.$ random variables $X_{1},\ldots,X_{n}$, we denote their summation by
$S_{n}\coloneqq\sum_{i=1}^{n}X_{i}$. When $n$ is even, we denote by
$S^{\mathsf{O}}_{n/2}\coloneqq\sum_{k=1}^{n/2}S_{2k-1}$ and
$S^{\mathsf{E}}_{n/2}\coloneqq\sum_{k=1}^{n/2}S_{2k}$ the summations of their
odd and even terms, respectively.

###  1.3 Assumptions

Throughout this paper, we assume that we can access a simulator $\cal S$. The
simulator can take any trajectory $y^{(0:d-1)}$ with $0\leq d\leq D$ as input,
and outputs $y^{(d)}$ which follows the distribution $\pi_{d}(\cdot\mid
y^{(0:d-1)})$. In particular, ${\mathcal{S}}$ can take $\varnothing$ as input
and simulates $y^{(0)}\sim\pi_{0}$. Calling ${\mathcal{S}}$ recursively for
$D+1$ times generates one complete sample path. This assumption enables us to
sample from any marginal or conditional distribution perfectly. This
assumption is also standard and is posed explicitly or implicitly in nearly
all the existing works concerning the estimation of nested expectations, see
Giles and Goda [2019], Goda et al. [2022], Zhou et al. [2022] for examples.

Fix a function $f:\mathbf{R}^{k+2}\rightarrow\mathbf{R}$, we say $f$ satisfies
the last-component bounded Lipschitz condition (LBL) if there exists
$L<\infty$ such that for all $x,z\in\mathbf{R}$:

$\displaystyle\sup_{y^{(0:k)}}\lvert f(y^{(0:k)},x)-f(y^{(0:k)},z)\rvert\leq
L\lvert x-z\rvert.$ (4)

We say $f$ satisfies the last-component bounded second-derivative condition
(LBS) if $f$ has continuous second-order derivative on its last component, and
there exists $C<\infty$ such that

$\displaystyle\sup_{y^{(0:k+1)}}\left|\partial^{2}_{k+1}f(y^{(0:k+1)})\right|<C,$
(5)

where $\partial^{2}_{k+1}f\coloneqq(\partial^{2}f)/(\partial{y^{(k+1)})^{2}}$
is the second-order partial derivative on the last component. These
assumptions (and their variants) are also posed in related works such as
Rainforth [2018], Blanchet and Glynn [2015], Giles [2018].

##  2 Algorithm, estimator, and theoretical results

Now we are ready to present our main results. As discussed in Section 1, we
will be focusing on designing a Monte Carlo estimator which is unbiased, has a
finite computational cost, and has finite variance or (2-$\delta$)-th moment
under different assumptions.

###  2.1 Preliminary analysis

One of the challenges in estimating the RNEs is the difficulty of estimating
$\gamma_{1}(y^{(0)})$. Users typically first estimate $\gamma_{1}(y^{(0)})$
and then use these estimators to estimate $\gamma_{0}$. For the time being, we
are temporarily adding the assumption that users can simulate unbiased
estimators $\hat{\gamma}_{1}(y^{(0)})$ of $\gamma_{1}(y^{(0)})$ for every
fixed $y^{(0)}$ with finite computational cost. This assumption will be
removed in Section 2.2. It easily holds when $D=1$, as users can repeatedly
simulate $y^{(1)}_{i}\sim\pi_{1}(\cdot\mid y^{(0)})$ and it follows from the
problem definition that each $g_{1}(y^{(0)},y^{(1)}_{i})$ is unbiased for
$\gamma_{1}(y^{(0)})$. In the general case of $D>1$, this assumption is far
from trivial, as $\gamma_{1}(y^{(0)})$ is itself a nested expectation with a
nesting depth of $D-1$. Nevertheless, as we will see in Section 2.2, this
assumption helps us to capture and reduce the intrinsic difficulty of the
problem and, therefore, will guide us to design the general algorithm.

With this extra assumption, constructing unbiased estimators of (1) is
equivalent to constructing unbiased estimators of
$g_{0}(y^{(0)},\gamma_{1}(y^{(0)}))$. Even with access to unbiased estimators
of $\gamma_{1}(y^{(0)})$, the intuitive plug-in estimator
$g_{0}\big{(}y^{(0)},\hat{\gamma}_{1}(y^{(0)})\big{)}$ is still biased, as in
general $\mathbf{E}[g_{0}\big{(}y^{(0)},\hat{\gamma}_{1}(y^{(0)})\big{)}\mid
y^{(0)}]\neq g_{0}(y^{(0)},\mathbf{E}[\hat{\gamma}_{1}(y^{(0)})\mid
y^{(0)}])$. To eliminate this bias, we use the rMLMC method Blanchet and Glynn
[2015], which is briefly reviewed below.

The rMLMC method uses the Law of Large Numbers (LLN) and rewrites $g_{0}$ as
the following telescoping summation.

$\displaystyle
g_{0}(y^{(0)},\gamma_{1}(y^{(0)}))=\mathbf{E}\left[g_{0}\left(y^{(0)},\lim_{k\rightarrow\infty}\frac{S_{k}}{k}\right)\bigg{|}~{}y^{(0)}\right]$
$\displaystyle=\sum_{n=1}^{\infty}\mathbf{E}\left[g_{0}\left(y^{(0)},\frac{S_{2^{n}}}{2^{n}}\right)\bigg{|}~{}y^{(0)}\right]$
$\displaystyle\quad\quad\quad\quad-\mathbf{E}\left[g_{0}\left(y^{(0)},\frac{S_{2^{n-1}}}{2^{n-1}}\right)\bigg{|}~{}y^{(0)}\right],$

where $S_{k}=\sum_{i=1}^{k}\hat{\gamma}_{1,i}(y^{(0)})$ is the summation of
i.i.d. copies of $\hat{\gamma}_{1}(y^{(0)})$. To unbiasedly estimate the
infinite sum, the rMLMC algorithm first samples $y^{(0)}\sim\pi_{0}$, then
samples a random $N\sim{\mathsf{Geo}}(r)$, finally generates $2^{N}$ unbiased
estimators $\\{\hat{\gamma}_{1,i}(y^{(0)})\\}_{i=1}^{2^{N}}$ of
$\gamma_{1}(y^{(0)})$ and estimates $\gamma_{0}$ by
$R_{0}\coloneqq\Delta_{N}/p_{r}(N)$, where $\Delta_{n}$ is defined as:

$\displaystyle\Delta_{n}\coloneqq
g_{0}\left(y^{(0)},\frac{S_{2^{n}}}{2^{n}}\right)-\frac{1}{2}\Bigg{[}$
$\displaystyle
g_{0}\bigg{(}y^{(0)},\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}\bigg{)}$
$\displaystyle+g_{0}\bigg{(}y^{(0)},\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}\bigg{)}\Bigg{]}$

for $n\geq 1$ and $\Delta_{0}\coloneqq
g_{0}(y^{(0)},\hat{\gamma}_{1,1}(y^{(0)}))$.

The next theorem justifies the theoretical properties of $R_{0}$:

######  Theorem 2.1.

With all the notations as above, suppose
$g_{0}:\mathbf{R}^{2}\rightarrow\mathbf{R}$ satisfies LBS condition defined in
(5), and $\lVert\hat{\gamma}_{1}(y^{(0)})\rVert_{\pi,m}<\infty$ for some
$m\geq 4$. Then for any $r\in(1/2,3/4)$, the estimator
$R_{0}\coloneqq\Delta_{N}/p_{r}(N)$ has expectation $\gamma_{0}$, finite
variance, and finite expected computational cost.

Theorem 2.1 will not be proved directly, as it is a special case of our
Theorem 2.2. For now, we use the following heuristic calculation to justify
the unbiasedness of $\hat{\gamma}_{0}$:

$\displaystyle\mathbf{E}[$ $\displaystyle
R_{0}\big{|}~{}y^{(0)}]=\mathbf{E}\left[\mathbf{E}\left[R_{0}\big{|}~{}N,y^{(0)}\right]\right]=\sum_{n=0}^{\infty}\mathbf{E}\left[\frac{\Delta_{n}}{p_{r}(n)}p_{r}(n)\big{|}~{}y^{(0)}\right]$
$\displaystyle=\sum_{n=0}^{\infty}\mathbf{E}\left[g_{0}\left(y^{(0)},\frac{S_{2^{n}}}{2^{n}}\right)\bigg{|}~{}y^{(0)}\right]-\mathbf{E}\left[g_{0}\left(y^{(0)},\frac{S_{2^{n-1}}}{2^{n-1}}\right)\bigg{|}~{}y^{(0)}\right]$
$\displaystyle=g_{0}(y^{(0)},\gamma_{1}(y^{(0)})).$

Therefore
$\mathbf{E}[R_{0}]=\mathbf{E}[g_{0}(y^{(0)},\gamma_{1}(y^{(0)}))]=\gamma_{0}$
by (1). More technical discussions such as the range of $r$, other possible
regularity conditions on $g_{0}$, and the moment guarantees of $\gamma_{0}$
will all be deferred after Theorem 2.2.

###  2.2 Recursive rMLMC algorithm for general $D$

Theorem 2.1 is useful to solve our original problem (without extra
assumptions) in two ways. First, Theorem 2.1 already solves the case where
$D=1$, as our extra assumption automatically holds. It states that if $g_{0}$
has a bounded second derivative on its last component, and
$g_{1}(y^{(0)},y^{(1)})$ has at least finite fourth moment under $\pi$, then
$R_{0}$ is unbiased, has finite variance, and finite expected computational
cost. More importantly, Theorem 2.1 tells us that the original problem of
estimating an RNE with a depth of $D$ can be solved if we can unbiasedly
estimate $\gamma_{1}(y^{(0)})$ for fixed $y^{(0)}$, which is another RNE with
a depth of $D-1$. Therefore, we have successfully reduced the number of
nestings by one. This observation motivates us to come up with an algorithm
for the general $D$ case, as explained below.

We first go one step further to illustrate the $D=2$ case. When $D=2$,
estimating $\gamma_{1}(y^{(0)})$ reduces to the case we have analyzed in
Section 2.1. To be precise, since $g_{2}(y^{(0:2)})$ is unbiased for
$\gamma_{2}(y^{(0:1)})$ if $y^{(2)}\sim\pi_{2}(\cdot\mid y^{(0:1)})$, one can
first sample $y^{(1)}\sim\pi_{1}(\cdot\mid y^{(0)})$, then simulate
$N\sim{\mathsf{Geo}}(r)$ and $2^{N}$ samples $\\{y^{(2)}_{i}\\}_{i=1}^{2^{N}}$
from $\pi_{2}(\cdot\mid y^{(0:1)})$. Let
$\hat{\gamma}_{2,i}(y^{(0:1)})\coloneqq g_{2}(y^{(0:1)},y^{(2)}_{i})$, our
estimator of $\gamma_{1}(y^{(0)})$ is then constructed in the same way as
Section 2.1, i.e., $R_{1}(y^{(0)}):=\Delta_{N}/p_{r}(N)$ with

$\displaystyle\Delta_{n}\coloneqq
g_{1}\left(y^{(0:1)},\frac{S_{2^{n}}}{2^{n}}\right)-\frac{1}{2}\Bigg{[}$
$\displaystyle
g_{1}\bigg{(}y^{(0:1)},\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}\bigg{)}$
$\displaystyle+g_{1}\bigg{(}y^{(0:1)},\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}\bigg{)}\Bigg{]},$

where $S_{2^{n}},S^{\mathsf{E}}_{2^{n-1}},S^{\mathsf{O}}_{2^{n-1}}$ are the
summation of every, even, and odd terms of
$\\{\hat{\gamma}_{2,i}(y^{(0:1)})\\}$, respectively. The same procedure of
simulating $R_{1}(y^{(0)})$ can be repeated independently. Therefore we can
sample another geometrically distributed random variable
$N^{\prime}\sim{\mathsf{Geo}}(r^{\prime})$, and generate
$R_{1,i}(y^{(0)}):=\Delta_{N^{\prime}}/p_{r}(N^{\prime})$ independently. Since
each $R_{1,i}(y^{(0)})$ is unbiased for $\gamma_{1}(y^{(0)})$, one can again
use the method described in Section 2.1 to form our final estimator for
$\gamma_{0}$. After checking $R_{1}(y^{(0)})$ satisfies the finite fourth-
moment assumption, Theorem 2.1 can be applied which implies our estimator is
unbiased, has finite variance and finite cost (for the $D=2$ case).

The general case works in the same way. A key observation is that, due to the
nested structure of the problem, Theorem 2.1 not only states that an unbiased
estimator of $\gamma_{0}$ can be constructed if one can unbiasedly estimate
$\gamma_{1}(y^{(0)})$ for every $y^{(0)}$, but also directly implies that an
unbiased estimator of $\gamma_{d}(y^{(0:d-1)})$ can be constructed if one can
unbiasedly estimate $\gamma_{d+1}(y^{(0:d)})$ for every $y^{(0:d)}$.
Therefore, we can estimate $\gamma_{0}$ in a backward, inductive manner.

To begin, we consider the deepest depth of the problem, fixing any
$y^{(0:D-1)}$. An unbiased estimator of $\gamma_{D}(y^{(0:D-1)})$ can be
directly constructed as $g_{D}(y^{(0:D-1)},y^{(D)})$, where
$y^{(D)}\sim\pi_{D}(\cdot\mid y^{(0:D-1)})$. For $0\leq d\leq D-1$, if we
assume that users can generate unbiased estimators of
$\gamma_{d+1}(y^{(0:d)})$ for every $y^{(0:d)}$, then we can obtain an
unbiased estimator of $\gamma_{d}(y^{(0:d-1)})$ by sampling one $y^{(d)}$,
generating $N_{d}\sim{\mathsf{Geo}}(r_{d})$ and $2^{N_{d}}$ unbiased
estimators of $\gamma_{d+1}(y^{(0:d)})$, and applying the method described in
Section 2.1. This process continues until we reach $d=0$, at which point we
have an unbiased estimator of $\gamma_{0}$. The parameters
$(r_{0},r_{1},\ldots,r_{D-1})$ will be carefully chosen and depend on the
regularity assumptions of $(g_{0},g_{1},\ldots,g_{D-1})$. These choices will
be discussed in more detail later.

Our algorithm is described in Algorithm 1. It is written as a recursive
algorithm, though it could also be equivalently written in an iterative form
with much more cumbersome notations. Algorithm 1 takes a depth index, a
trajectory, a simulator, and parameters for the geometric distribution as
inputs, and outputs an unbiased estimator of $\gamma_{d}(H)$. In particular,
with inputs {depth $=0$, trajectory = $\varnothing$, parameters =
$(r_{0},r_{1},\ldots,r_{D-1})$}, it outputs $\mathsf{READ}$ – an unbiased
estimator of the RNE defined in (1). The logic of Algorithm 1 is precisely the
same as we just discussed. To estimate $\gamma_{d}(y^{(0:d-1)})$, the
algorithm first checks the value of $d$. When $d=D$, the problem becomes
straightforward. When $d<D$, the algorithm samples $y^{(d)}$, appends
$y^{(d)}$ to the trajectory, samples $N_{d}$, and calls itself $2^{N_{d}}$
times with depth $d+1$ and new trajectory $\\{y^{(0:d)}\\}$ to get $2^{N_{d}}$
unbiased estimators of $\gamma_{d+1}(y^{(0:d)})$. Finally, we split these
$2^{N_{d}}$ estimators into even and odd terms and apply the method described
in Section 2.1. The algorithm is guaranteed to stop, as the depth will
eventually reach the deepest depth $D$.

Algorithm 1 A recursive rMLMC algorithm for RNEs

Input: Depth index $d\in\\{0,...,D\\}$. Trajectory history
$H=\\{y^{0},...,y^{d-1}\\}$ or $\varnothing$. A simulator $\mathcal{S}$.
Parameters $r_{d},...,r_{D-1}$ determined by conditions on
$\\{g_{d}\\}_{d=0}^{D-1}$ (see Theorem 2.2, 2.4).

Output: An unbiased estimator of $\gamma_{d}(H)$

if $d=D$ then

Sample one $y^{(D)}\sim\pi_{D}\left(\cdot\mid H\right)$;

Return: $R_{D}\coloneqq g_{D}\left(y^{(0:D)}\right).$

else

Sample $y^{(d)}\sim\pi_{d}\left(\cdot\mid H\right)$;

Update the trajectory $H\leftarrow H\cup\left\\{y^{(d)}\right\\}$;

Sample $N_{d}\sim{\mathsf{Geo}}(r_{d})$;

Call Algorithm 2 for $2^{N_{d}}$ times with inputs
$\\{d+1;H;\mathcal{S};r_{d+1},...,r_{D-1}\\}$, and label the observations as
$R_{d+1}(y^{(0:d)})(1),...,R_{d+1}(y^{(0:d)})\left(2^{N_{d}}\right)$;

Calculate
$S_{2^{N_{d}}},S^{\mathsf{E}}_{2^{N_{d}-1}},S^{\mathsf{O}}_{2^{N_{d}-1}}$
defined in Section 1.2;

Calculate $\left(\text{note }\Delta_{0}\coloneqq
g_{d}\left(y^{(0:d)},R_{d+1}(y^{(0:d)})(1)\right)\right)$:
$\displaystyle\Delta_{N_{d}}$
$\displaystyle=g_{d}\left(y^{(0:d)},\frac{S_{2^{N_{d}}}}{2^{N_{d}}}\right)-$
$\displaystyle\frac{1}{2}\left[g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{O}}_{2^{{N_{d}}-1}}}{2^{{N_{d}}-1}}\right)+g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{E}}_{2^{{N_{d}}-1}}}{2^{{N_{d}}-1}}\right)\right];$

Return: $R_{d}\coloneqq\Delta_{N_{d}}/p_{r_{d}}(N_{d})$.

end if

###  2.3 Theoretical guarantees

We now discuss the computational costs of Algorithm 1 and the statistical
properties of $\mathsf{READ}$. Our theoretical results depend on the
smoothness conditions of $\\{g_{d}\\}_{d=0}^{D-1}$, so we will consider two
cases where $\\{g_{d}\\}_{d=0}^{D-1}$ satisfies the LBS and LBL assumptions
separately.

####  2.3.1 The LBS case

The following theorem shows, under the LBS assumption, the computational cost
and the variance of $\mathsf{READ}$ can be controlled simultaneously.

######  Theorem 2.2.

Suppose for every $d\in\\{0,1,\ldots,D-1\\}$, the function $g_{d}$ satisfies
the LBS assumption defined in (5), and $r_{d}\coloneqq 1-2^{-k_{d}}$ satisfies
$k_{d}\in\left(1,\frac{2^{d+1}}{2^{d+1}-1}\right)$. Moreover, suppose $\lVert
g_{D}(y^{(0:D)})\rVert_{\pi,2^{D+1}}<\infty$. Then for every $0\leq d\leq D$,
the output $R_{d}(y^{(0:d-1)})$ of Algorithm 1 with inputs {depth = $d$,
trajectory = $y^{(0:d-1)}$, ${\mathcal{S}}$, parameters
$(r_{d},\ldots,r_{D-1})$} has the following properties:

  * •

For $\pi$-almost surely every fixed $y^{(0:d-1)}$,

$\mathbf{E}\left[R_{d}(y^{(0:d-1)})\mid
y^{(0:d-1)}\right]=\gamma_{d}(y^{(0:d-1)}).$

  * •

The expected computational cost of $R_{d}$ is finite.

  * •

The output has finite $2^{d+1}$-th moment, i.e.,

$\mathbf{E}_{\pi}\left[|R_{d}(y^{(0:d-1)})|^{2^{d+1}}\right]<\infty~{}~{}\text{for}~{}~{}0\leq
d\leq D.$

Theorem 2.2 states for $\pi$-almost surely every $y^{(0:d-1)}$, the
expectation of the output $R_{d}$ conditioning on the input is unbiased for
$\gamma_{d}(y^{(0:d-1)})$. The computational cost has a finite expectation,
and the output has a finite $2^{d+1}$-th moment 111Readers should notice that
the expectation of $R_{d}(y^{(0:d-1)})$ is calculated under the conditional
distribution $\pi_{d:D}(\cdot\mid y^{(0:d-1)})$. The computational cost and
the $2^{d+1}$-th moment are calculated under the joint distribution $\pi$.
When the input depth $=0$, these two underlying distributions coincide.. The
detailed proof of Theorem 2.2 will be provided in the Appendix. Here, we
highlight two special cases. First, Theorem 2.2 shows that $\mathsf{READ}$,
the output $R_{0}$ of Algorithm 1 when given input {depth = $0$, trajectory =
$\varnothing$, ${\mathcal{S}}$, parameters = $(r_{0},\ldots,r_{D-1})$}, has
the desired properties. Specifically, it is an unbiased estimator for
$\gamma_{0}$ with finite expected computational cost and finite variance.
Second, Theorem 2.2 includes Theorem 2.1 as a special case when $D=1$.

Let $R_{0,1},R_{0,2},\ldots$ be the i.i.d. outcomes by repeatedly implementing
Algorithm 1. Since each $R_{0,i}$ is unbiased and has a finite variance, the
standard Central Limit Theorem (CLT) implies that
$\sqrt{n}(\sum_{i=1}^{n}R_{0,i}/n-\gamma_{0})\rightarrow\mathbf{N}(0,1)$ in
distribution. This means that the estimator $\sum_{i=1}^{n}R_{0,i}/n$
converges to $\gamma_{0}$ at a rate of $n^{-1/2}$ in rMSE (or equivalently,
$n^{-1}$ in MSE), which compares quite favorably with the rates obtained by
NMC estimators in Rainforth [2018]. This rate is optimal in the sense that it
matches the minimax lower bound over all Monte Carlo methods (Theorem 2.1 of
Heinrich and Sindambiwe [1999]). The next corollary shows that, by repeatedly
implementing Algorithm 1, users can easily obtain an unbiased estimator for
$\gamma_{0}$ with at most $\varepsilon$-rMSE within $O(\varepsilon^{-2})$
computational cost.

######  Corollary 2.3.

With all the assumptions the same as Theorem 2.2, for any $\varepsilon>0$, we
can construct an estimator $R$ with expected computational cost
${\mathcal{O}}(\varepsilon^{-2})$ such that the rMSE
$\sqrt{\mathbf{E}[(R-\gamma_{0})^{2}]}$ is at most $\epsilon$.

###### Proof of Corollary 2.3.

Calling Algorithm 1 independently for $n$ times with {depth = $0$, trajectory
= $\varnothing$, ${\mathcal{S}}$, parameters = $(r_{0},\ldots,r_{D-1})$} yield
i.i.d. unbiased estimators $R_{0,1},...,R_{0,n}$ for $\gamma_{0}$. Let our
estimator be $R\coloneqq\frac{1}{n}\sum_{i=1}^{n}R_{0,i}$. Then,

$\mathbf{E}[(R-\gamma_{0})^{2}]=\mathbf{E}\left[\left(\frac{1}{n}\sum_{i=1}^{n}R_{0,i}-\gamma_{0}\right)^{2}\right]=\frac{1}{n}\text{Var}(R_{0}).$

Thus noting that $\text{Var}(R_{0})<\infty$ by Theorem 2.2, taking
$n=\text{Var}(R_{0})/\varepsilon^{2}$ samples ensures $R$ has up to
$\varepsilon$-rMSE. Finally, let $C\coloneqq C(D)<\infty$ be the expected
computational cost for one call of Algorithm 1 . The expected computational
cost for constructing $R$ is then
$C\cdot\text{Var}(R_{0})/\varepsilon^{2}=\Theta(\varepsilon^{-2})$. ∎

The $\varepsilon$-rMSE of $R$ can be easily translated to other performance
metrics via standard inequalities. For example, for any $\delta$, Markov’s
inequality implies the absolute error $\lvert R-\gamma_{0}\rvert$ is less than
$\varepsilon/\sqrt{\delta}$ with probability at least $1-\delta$.

Next, we discuss the assumptions and proof strategies of Theorem 2.2. We
require the first $D$ functions $\\{g_{d}\\}_{d=0}^{D-1}$ all satisfy the LBS
condition, and the final function $g_{D}$ has finite $2^{D+1}$-th moment under
$\pi$. The LBS assumption also appears in the work of NMC estimators (see the
second part of Theorem 3 in Rainforth et al. [2018]). The moment assumption of
$g_{D}$ is not required in Rainforth et al. [2018]. Nevertheless, it is a mild
assumption that holds in most practical applications. It covers all the cases
where $g_{D}$ is bounded or has a moment generating function (including the
uniform, Gaussian, Poisson, or exponential distributions), which implies
$\mathbf{E}[\lvert g_{D}\rvert^{k}]<\infty$ for every $k$. As we will see in
our proofs, these assumptions help us to establish the moment guarantee of
Theorem 2.2 in a backward inductive way. For example, the $2^{D+1}$-th moment
assumption on $g_{D}$ and the LBS assumption on $g_{D-1}$ implies $R_{D-1}$
has finite $2^{D}$-th moment. More generally, the finiteness of the
$2^{d+1}$-th moment of $R_{d}$ follows from the LBS assumption on $g_{d}$ and
the $2^{d+2}$-th moment of $R_{d+1}$ (which is the conclusion of the previous
inductive step). Eventually, we conclude $R_{0}$ has a finite variance.
Finally, we want to emphasize our moment assumption on $g_{D}$ is not
‘trajectory-dependent’. We require $g_{D}$ has finite $2^{D+1}$-th moment
under the joint distribution $\pi$ of $y^{(0:D)}$, which is much weaker than
$g_{D}$ has a uniformly bounded finite $2^{D+1}$-th moment under
$\pi_{D}(\cdot\mid y^{(0:D-1)})$ for every fixed trajectory $y^{(0:D-1)}$.

Finally, the parameters $\\{r_{d}\\}_{d=0}^{D=1}$ reflect the trade-off
between the variance and computation cost. Since $2^{N_{d}}$ calls are
required for each $d$, standard calculation shows that
$\mathbf{E}[2^{N_{d}}]=r_{d}/(2r_{d}-1)$ when $r_{d}>0.5$, and $+\infty$ if
$r_{d}\leq 0.5$. Therefore, every $r_{d}$ has to be strictly greater than
$0.5$ to ensure a finite expected computational cost. Meanwhile, we cannot
guarantee finite variance or unbiasedness of $\mathsf{READ}$ when $r_{d}$
becomes too large. Our range for $r_{d}$ in Theorem 2.2 follows from a careful
calculation in our proof to ensure unbiasedness, finite computational cost,
and variance simultaneously.

####  2.3.2 The LBL case

The assumptions in Theorem 2.2 guarantee that $\mathsf{READ}$ enjoys an
optimal convergence rate and computational cost. However, the second-order
derivative assumption also rules out many functions of practical interest,
such as $\max$ and $\min$. In this section, we study the theoretical
properties of Algorithm 1 and $\mathsf{READ}$ under weaker smoothness and
moment assumptions. Our result is summarized below:

######  Theorem 2.4.

Fix any $0<\delta<1/2$. Suppose for every $d\in\\{0,1,\ldots,D-1\\}$, the
function $g_{d}$ satisfies the LBL assumption defined in (4), and
$r_{d}\coloneqq 1-2^{-k_{d}}$ satisfies

$k_{d}\in\left(1,\left(\frac{2^{d+2}-3\delta}{2^{d+3}-3\delta}\right)\left(\frac{2^{d+1}-\delta}{2^{d}-\delta}\right)\right).$

Moreover, suppose $\lVert g_{D}(y^{(0:D)})\rVert_{\pi,2}<\infty$. Then for
every $0\leq d\leq D$, the output $R_{d}(y^{(0:d-1)})$ of Algorithm 1 with
inputs {depth = $d$, trajectory = $y^{(0:d-1)}$, ${\mathcal{S}}$, parameters
$(r_{d},\ldots,r_{D-1})$} has the following properties:

  * •

For almost surely every fixed $y^{(0:d-1)}$,

$\mathbf{E}\left[R_{d}(y^{(0:d-1)})\mid
y^{(0:d-1)}\right]=\gamma_{d}(y^{(0:d-1)}).$

  * •

The expected computational cost of $R_{d}$ is finite.

  * •

The output has finite $(2-\delta/2^{d})$-th moment, i.e.,

$\mathbf{E}_{\pi}\left[|R_{d}(y^{(0:d-1)})|^{(2-\delta/2^{d})}\right]<\infty~{}~{}\text{for}~{}~{}0\leq
d\leq D.$

Comparing Theorem 2.2, which requires the LBS assumption for
$\\{g_{d}\\}_{d=0}^{D-1}$ and finite $2^{D+1}$-th moment for $g_{D}$, with
Theorem 2.4, which only requires the LBL assumption for
$\\{g_{d}\\}_{d=1}^{D-1}$ and finite second moment for $g_{D}$, it is clear
that Theorem 2 has more general assumptions. However, it does not guarantee
that $\mathsf{READ}$ has a finite variance. Nevertheless, it remains unbiased
and has a finite expected computational cost. To minimize the loss of moment
guarantees, one can choose suitable parameters such that $\mathsf{READ}$ has
finite $(2-\delta)$-th moment for any small $\delta$.

Again, let $R_{0,1},R_{0,2},\ldots,$ be the i.i.d. outcomes by repeatedly
implementing Algorithm 1. There are more technical challenges when analyzing
the convergence rate of $\sum_{i=1}^{n}R_{0,i}/n-\gamma_{0}$, as the CLT
cannot be applied. Instead, we use the Marcinkiewicz-Zygmund generalized law
of large numbers (see Theorem A.4 in Appendix A), which shows
$n^{-1}\mathbf{E}[\lvert\sum_{i=1}^{n}X_{i}\rvert^{p}]\rightarrow 0$ if
$\\{X_{i}\\}_{i=1}^{n}$ are i.i.d., centered random variables with finite
$p$-th moment for $p\in[1,2)$. Our result is the following:

######  Corollary 2.5.

With all the assumptions the same as Theorem 2.4, let
$R_{0,1},R_{0,2},\ldots,$ be the i.i.d. outcomes by repeatedly implementing
Algorithm 1, we have:

  * •

$\mathbf{E}[\lvert\sum_{i=1}^{n}R_{0,i}/n-\gamma_{0}\rvert]=o(n^{-1/(2(1+\delta))})$.

  * •

We can construct an estimator $R$ with expected computational cost
${\mathcal{O}}(\varepsilon^{-2(1+\delta)})$ such that the mean absolute error
$\mathbf{E}[\lvert R-\gamma_{0}\rvert]<\varepsilon$.

###### Proof of Corollary 2.5.

Applying Theorem A.4 with $p=2-\delta,X_{i}=R_{0,i}-\gamma_{0}$ and Jensen’s
inequality, we have:

$\displaystyle\mathbf{E}\left[\left\lvert\sum_{i=1}^{n}R_{0,i}/n-\gamma_{0}\right\rvert\right]=n^{-1}\mathbf{E}\left[\left\lvert\sum_{i=1}^{n}X_{i}\right\rvert\right]$
$\displaystyle\leq
n^{-1}\left(\mathbf{E}\left[\left\lvert\sum_{i=1}^{n}X_{i}\right\rvert^{p}\right]\right)^{1/p}=o(n^{-1+\frac{1}{p}})=o(n^{\frac{-1}{2(1+\delta)}}),$

which proves the first part. The last step follows from
$1-1/(2-\delta)>1/(2+2\delta)$ for $\delta\in(0,1/2)$. Setting
$n=\Omega(\varepsilon^{-2(1+\delta)})$ and the second part immediately
follows. ∎

Although we are not able to recover the optimal $n^{-1/2}$ convergence rate
under this more general regime, our convergence rate is still near-optimal as
it can be as close to $n^{-1/2}$ as we want. Still, the convergence rate does
not depend on $D$, and, although we replace the MSE by MAE due to the moment
constraint, one can still use Markov’s inequality to show the absolute error
$\lvert R-\gamma_{0}\rvert$ is less than $\varepsilon/\delta$ with probability
at least $1-\delta$.

As the $\max$ function satisfies the LBL assumption, our results here include
the optimal stopping problem as a special case. Our results complement the
work of Zhou et al. [2022], where the authors use rMLMC to design an estimator
with ${\mathcal{O}}(\varepsilon^{-2})$ computational cost under stronger
assumptions (see their Assumption 4). We have a slightly worse cost of
${\mathcal{O}}(\varepsilon^{-2(1+\delta)})$ but under more general
assumptions.

##  3 Numerical experiments

We consider the following simple example with known ground-truth. Suppose the
process $(y^{(0)},y^{(1)},y^{(2)})$ satisfies
$y^{(0)}\sim\mathbf{N}(\pi/2,1),y^{(1)}\sim\mathbf{N}(y^{(0)},1),y^{(2)}\sim\mathbf{N}(y^{(1)},1)$.
Define
$g_{0}(y^{(0)},z)\coloneqq\sin(y^{(0)}+z),g_{1}(y^{(0:1)},z)\coloneqq\sin(y^{(1)}-z)$,
and $g_{2}(y^{(0:2)})\coloneqq y^{(2)}$. The target quantity $\gamma_{0}$
defined (1) is a nested expectation with $D=2$. One can use the formula
$\mathbf{E}_{Z\sim\mathbf{N}(\mu,\sigma^{2})}[\sin(Z)]=\sin(\mu)\exp(-\sigma^{2}/2)$
to analytically calculate $\gamma_{0}=\exp(-1/2)\approx 0.6065$. Now we
compare our $\mathsf{READ}$ estimator with the NMC estimators in Rainforth et
al. [2018].

For the NMC estimator, users first specify $N_{0},N_{1},N_{2}$. Then we sample
$N_{0}$ copies of $y^{(0)}$, $N_{1}$ copies of $y^{(1)}$ for each fixed
$y^{(0)}$, and $N_{2}$ copies of $y^{(2)}$ for each fixed $y^{(0:1)}$, and use
these samples to form the NMC estimator, details are explained in Appendix D.
Following Rainforth [2018], we consider two ways of allocating
$(N_{0},N_{1},N_{2})$. The first estimator NMC1 is to choose
$N_{0}=N_{1}=N_{2}$, the second NMC2 is to choose $N_{0}=N_{1}^{2}=N_{2}^{2}$.
For $\mathsf{READ}$, since all assumptions in Theorem 2.2 are satisfied,
therefore when $r_{0}\in(1/2,3/4)$ and $r_{1}\in(1/2,1-2^{-4/3})$, the
$\mathsf{READ}$ estimator generated by Algorithm 1 is unbiased and of finite
variance. Since the computational cost gets lower when each $r_{i}$ gets
larger, we choose $r_{0}=0.74$ and $r_{1}=0.6$ (close to the upper-end of
their respective ranges above) to facilitate the computational efficiency.
Therefore, implementing Algorithm 1 once has an expected sample
size/computational cost
$\left(r_{1}/(2r_{1}-1)\right)\left(r_{2}/(2r_{2}-1)\right)\approx 4.625$.

Our comparison result is summarized in Figure 1. Since the NMC methods and
$\mathsf{READ}$ have different ways of generating estimators. To make a fair
comparison, we compare the estimation errors with the total sample cost used
by these three estimators. For NMC1 and NMC2, the total sample cost is
$n=N_{0}N_{1}N_{2}$. For $\mathsf{READ}$, the total sample cost is random,
therefore we use its expected value, which equals $4.625\times$ Number of
repetitions of Algorithm 1. The slopes of the blue, red, and green lines,
which correspond to the empirical convergence rate of $\mathsf{READ}$, NMC1,
NMC2, equals $-0.97,-0.35,-0.47$, respectively. They match well with the
theoretical predictions $n^{-1}$ in Corollary 2.3 for $\mathsf{READ}$,
$n^{-1/3}$ for NMC1, and $n^{-1/2}$ for NMC2 in Rainforth et al. [2018]. It is
clear from Figure 1 that READ has a significant advantage over NMC estimators,
with both faster convergence rate and orders of magnitude lower errors.

Figure 1: The comparison on the empirical MSEs of estimating the RNE among
$\mathsf{READ}$ (blue), NMC1 (red), and NMC2 (green). All the logarithms are
of the base $10$. Each method’s empirical errors are calculated based on $20$
independent repetitions.

We also repeatedly call Algorithm 1 for $10^{6}$ times and plot the running
averages of our estimates in Figure 2. Our estimator becomes more accurate
when we increase the number of repetitions. For each
$k\in(1,2,\ldots,10^{6})$, we also calculate the standard deviation (SD) of
the first $k$ repetitions and use Mean $\pm 1.96$ SD to form the $95\%$
confidence interval. It is also clear from Figure 2 that our confidence
intervals always include the ground-truth, suggesting the high accuracy of our
method. In contrast, constructing confidence intervals of NMC estimators are
much more time-consuming.

Some additional statistics and an extra experiment are provided in Appendix E
and F.

Figure 2: The trace plot (solid blue curve) of the running averages of
$\mathsf{READ}$. The blue dotted curves are the $95\%$ confidence intervals.
The red dashed line is the ground truth $\exp(-1/2)$.

##  4 Further discussions

Here we provide some remarks for practical implementation and discuss some
potential generalizations. The users need to specify the parameters
$\\{r_{d}\\}_{d=0}^{D-1}$ when implementing Algorithm 1. Larger values of
$r_{i}$ lead to a shorter time for each implementation but potentially larger
variance. When some $r_{i}$ is not chosen according to Theorem 2.2 or 2.4, the
algorithm can still be implemented, but the variance may be infinite. The
trade-off between the values of $\\{r_{d}\\}$ and the fluctuations of the
resulting estimator is problem-specific. In practice, knowing how many
repetitions are sufficient is important to provide an accurate estimator. One
possible way is to bound certain moments of $\mathsf{READ}$ and use Corollary
2.3 or 2.5 to choose a sufficiently large $n$. But this bound can be problem-
specific and very conservative. Instead, we follow Glynn and Rhee [2014] and
suggest the following adaptive stopping rule: users first specify a precision-
level $\varepsilon$ and a small $\delta\%$. When repeatedly implementing
Algorithm 1, users calculate the empirical $(1-\delta\%)$ confidence interval
$[L_{\delta}(k),U_{\delta}(k)]$ for first $k$ repetitions in the same way as
Section 3 for every $k$. Users can stop when the width of the confidence
interval is less than $2\varepsilon$. The validity of this stopping rule is
proven in Glynn and Whitt [1992].

There are two directions for further extensions. First, in our paper, each
$y^{(i)}$ is assumed to be univariate for concreteness. Extending this
assumption to multivariate $y^{(i)}$ should be straightforward and the
computational cost of Algorithm 1 scales linearly with the dimensionality of
$y^{(i)}$, which may be another appealing property of our method. Second, in
this paper, we only consider the ‘fix $D$’ regime and construct estimators
with optimality guarantees. However, the cost of Algorithm 1 scales
exponentially with $D$. Therefore, although our algorithm is more efficient
than the NMC estimator for every fixed $D$, both methods are not practical
when $D$ becomes too large. Indeed, the poor scaling with $D$ is a common
issue in related literature such as Glasserman and Yu [2004], Zanger [2013]
and seems unavoidable. An interesting direction would be to construct
modifications of Algorithm 1 under extra practical assumptions for large or
infinite $D$. For example, if we know that the ‘influence’ of $\gamma_{d}$ on
$\gamma_{0}$ decays exponentially or double-exponentially with $d$, it is then
sufficient to truncate the depth to $\tilde{D}\coloneqq\log(1/\varepsilon)$ or
$\sqrt{\log(1/\varepsilon)}$. We hope to report progress in the future.

## References

  * Biswas et al. [2019] N. Biswas, P. E. Jacob, and P. Vanetti.  Estimating convergence of Markov chains with L-lag couplings.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Blanchet and Glynn [2015] J. H. Blanchet and P. W. Glynn.  Unbiased Monte Carlo for optimization and functions of expectations via multi-level randomization.  _2015 Winter Simulation Conference (WSC)_ , pages 3656–3667, 2015\.
  * Broadie et al. [2011] M. Broadie, Y. Du, and C. C. Moallemi.  Efficient risk estimation via nested sequential simulation.  _Management Science_ , 57(6):1172–1194, 2011\.
  * Dauchet et al. [2018] J. Dauchet, J.-J. Bezian, S. Blanco, C. Caliot, J. Charon, C. Coustet, M. El Hafi, V. Eymet, O. Farges, V. Forest, et al.  Addressing nonlinearities in Monte Carlo.  _Scientific reports_ , 8(1):1–11, 2018.
  * Giles [2008] M. B. Giles.  Multilevel monte carlo path simulation.  _Operations research_ , 56(3):607–617, 2008.
  * Giles [2018] M. B. Giles.  MLMC for nested expectations.  In _Contemporary Computational Mathematics-A Celebration of the 80th Birthday of Ian Sloan_ , pages 425–442. Springer, 2018.
  * Giles and Goda [2019] M. B. Giles and T. Goda.  Decision-making under uncertainty: using MLMC for efficient estimation of EVPPI.  _Statistics and computing_ , 29(4):739–751, 2019\.
  * Giles and Haji-Ali [2019] M. B. Giles and A.-L. Haji-Ali.  Multilevel nested simulation for efficient risk estimation.  _SIAM/ASA Journal on Uncertainty Quantification_ , 7(2):497–525, 2019.
  * Giles and Haji-Ali [2022] M. B. Giles and A.-L. Haji-Ali.  Multilevel path branching for digital options.  _arXiv preprint arXiv:2209.03017_ , 2022.
  * Glasserman and Yu [2004] P. Glasserman and B. Yu.  Number of paths versus number of basis functions in American option pricing.  _The Annals of Applied Probability_ , 14(4):2090–2119, 2004.
  * Glynn and Rhee [2014] P. W. Glynn and C.-h. Rhee.  Exact estimation for Markov chain equilibrium expectations.  _Journal of Applied Probability_ , 51(A):377–389, 2014.
  * Glynn and Whitt [1992] P. W. Glynn and W. Whitt.  The asymptotic validity of sequential stopping rules for stochastic simulations.  _The Annals of Applied Probability_ , 2(1):180–198, 1992.
  * Goda et al. [2022] T. Goda, T. Hironaka, W. Kitade, and A. Foster.  Unbiased MLMC stochastic gradient-based optimization of bayesian experimental designs.  _SIAM Journal on Scientific Computing_ , 44(1):A286–A311, 2022.
  * Gordy and Juneja [2010] M. B. Gordy and S. Juneja.  Nested simulation in portfolio risk measurement.  _Management Science_ , 56(10):1833–1848, 2010\.
  * Gut [2005] A. Gut.  _Probability: A Graduate Course_.  Springer, 2005.
  * Heinrich [2001] S. Heinrich.  Multilevel Monte Carlo methods.  In _International Conference on Large-Scale Scientific Computing_ , pages 58–67. Springer, 2001.
  * Heinrich and Sindambiwe [1999] S. Heinrich and E. Sindambiwe.  Monte Carlo complexity of parametric integration.  _Journal of Complexity_ , 15(3):317–341, 1999\.
  * Jacob et al. [2020] P. E. Jacob, J. O’Leary, and Y. F. Atchadé.  Unbiased Markov chain Monte Carlo methods with couplings.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82(3):543–600, 2020.
  * Kahale [2022] N. Kahale.  Unbiased time-average estimators for Markov chains.  _arXiv preprint arXiv:2209.09581_ , 2022.
  * McLeish [2011] D. McLeish.  A general method for debiasing a Monte Carlo estimator.  _Monte Carlo methods and applications_ , 17(4):301–315, 2011.
  * Rainforth [2018] T. Rainforth.  Nesting probabilistic programs.  In _Conference on Uncertainty in Artificial Intelligence_ , pages 249–258, 2018.
  * Rainforth et al. [2018] T. Rainforth, R. Cornish, H. Yang, A. Warrington, and F. Wood.  On nesting Monte Carlo estimators.  In _International Conference on Machine Learning_ , pages 4267–4276. PMLR, 2018.
  * Rhee and Glynn [2015] C.-H. Rhee and P. W. Glynn.  Unbiased estimation with square root convergence for SDE models.  _Operations Research_ , 63(5):1026–1043, 2015\.
  * Vihola [2018] M. Vihola.  Unbiased estimators and multilevel Monte Carlo.  _Operations Research_ , 66(2):448–462, 2018.
  * Wang and Wang [2022] G. Wang and T. Wang.  Unbiased Multilevel Monte Carlo methods for intractable distributions: MLMC meets MCMC.  _arXiv preprint arXiv:2204.04808_ , 2022.
  * Wang et al. [2021] G. Wang, J. O’Leary, and P. Jacob.  Maximal couplings of the metropolis-hastings algorithm.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1225–1233. PMLR, 2021.
  * Zanger [2013] D. Z. Zanger.  Quantitative error estimates for a least-squares Monte Carlo algorithm for American option pricing.  _Finance and Stochastics_ , 17(3):503–534, 2013\.
  * Zhou et al. [2022] Z. Zhou, G. Wang, J. H. Blanchet, and P. W. Glynn.  Unbiased optimal stopping via the MUSE.  _Stochastic Processes and their Applications_ , 2022.

##  Appendix A Auxiliary results

The following theorem is from pg. 150 [Gut, 2005].

######  Theorem A.1.

Let $p\geq 1$. Suppose that $X_{1},X_{2},...X_{n}$ are independent random
variables. with mean 0 such that $\mathbf{E}|X_{k}|^{p}<\infty$ for all $k$,
and let $S_{n}\coloneqq\sum_{i=1}^{n}X_{i}$ denote the partial sums. Then
there exist constants $A_{p}^{*},B_{p}^{*}$ depending only on $p$ such that

$A_{p}^{*}\mathbf{E}\left(\sum_{k=1}^{n}X_{k}^{2}\right)^{p/2}\leq\mathbf{E}|S_{n}|^{p}\leq
B_{p}^{*}\mathbf{E}\left(\sum_{k=1}^{n}X_{k}^{2}\right)^{p/2}.$

The following corollary to the above theorem is from pg. 151 [Gut, 2005],
Corollary 8.2.

######  Corollary A.2.

Let $p\geq 1$. Suppose that $X,X_{1},X_{2}...X_{n}$ are i.i.d. random
variables. with mean 0 such that $\mathbf{E}|X|^{p}<\infty$, and let
$S_{n}\coloneqq\sum_{i=1}^{n}X_{i}$ denote the partial sums. Then there exists
a constant $B_{p}$ depending only on $p$, such that

$\mathbf{E}|S_{n}|^{p}\leq\begin{cases}B_{p}n^{p/2}\mathbf{E}|X|^{p},&p>2\\\
B_{p}n\mathbf{E}|X|^{p},&1\leq p\leq 2.\end{cases}$

The following lemma is instrumental in the proofs for the theoretical
guarantees of our algorithm under both the LBS and LBL assumptions.

######  Lemma A.3.

Let $(Z_{1},Z_{2})$ be a 2-stage stochastic process and there exists $p\geq
1$, such that $\mathbf{E}[|Z_{2}|^{p}]<\infty$. Conditioning on $Z_{1}$,
sample $i.i.d.$ $Z_{2}(1),...,Z_{2}(n)$. Then,

$\mathbf{E}\left[\left|\frac{1}{n}\sum_{i=1}^{n}Z_{2}(i)-\mathbf{E}[Z_{2}\mid
Z_{1}]\right|^{p}\right]\leq\begin{cases}B_{p}^{\prime}\frac{\mathbf{E}[|Z_{2}|^{p}]}{n^{p/2}}&p>2\\\
B_{p}^{\prime}\frac{\mathbf{E}[|Z_{2}|^{p}]}{n^{p-1}}&1\leq p\leq
2\end{cases}$

###### Proof.

Let $p>2$. For arbitrary fixed $Z_{1}=z_{1}$, define $\bar{Z}_{2}(i)\coloneqq
Z_{2}(i)-\mathbf{E}[Z_{2}(i)\mid Z_{1}=z_{1}]$, and apply Corollary A.2 on the
$i.i.d.$ mean $0$ random variables $\bar{Z}_{2}(i)$ under the probability
distribution $\pi(\cdot\mid Z_{1}=z_{1})$, we have

$\displaystyle\mathbf{E}\left[\left|\frac{1}{n}\sum_{i=1}^{n}Z_{2}(i)-\mathbf{E}[Z_{2}\mid
Z_{1}]\right|^{p}\right]$
$\displaystyle=\int_{\Omega}\mathbf{E}\left[\left|\frac{1}{n}\sum_{i=1}^{n}Z_{2}(i)-\mathbf{E}[Z_{2}\mid
Z_{1}]\right|^{p}\mid Z_{1}=z_{1}\right]\pi_{1}(\mathrm{d}z_{1})$
$\displaystyle=\frac{1}{n^{p}}\mathbf{E}[\left|\sum_{i=1}^{n}\bar{Z}_{2}(i)\right|^{p}\mid
Z_{1}=z_{1}]\pi_{1}(\mathrm{d}z_{1})$
$\displaystyle\leq\frac{B_{p}}{n^{p/2}}\mathbf{E}[\lvert\bar{Z}_{2}(1)\rvert^{p}\mid
Z_{1}=z_{1}]\pi_{1}(\mathrm{d}z_{1})$
$\displaystyle\leq\frac{B_{p}^{\prime}}{n^{p/2}}\mathbf{E}[\lvert
Z_{2}(1)\rvert^{p}\mid Z_{1}=z_{1}]\pi_{1}(\mathrm{d}z_{1})$
$\displaystyle=\frac{B_{p}^{\prime}}{n^{p/2}}\mathbf{E}[\lvert
Z_{2}(1)\rvert^{p}\rvert].$

The second inequality follows from the inequality $(a+b)^{p}\leq
2^{p-1}(\lvert a\rvert^{p}+\lvert b\rvert^{p})$ and the monotonicity of a
random variable’s $L^{p}$ norm :

$\displaystyle\mathbf{E}[|X-\mathbf{E}[X]|^{p}]\leq
2^{p-1}(\mathbf{E}[|X|^{p}]+\lvert\mathbf{E}[X]\rvert^{p})\leq
2^{p}\mathbf{E}[\lvert X\rvert^{p}]$

For the case $1\leq p\leq 2$, the calculation is identical to the above,
except replace the $B_{p}^{\prime}/n^{p/2}$ with $B_{p}^{\prime}/n^{p-1}$ from
Corollary A.2. ∎

The following theorem is the Marcinkiewicz-Zygmund law of large numbers from
pg. 311 [Gut, 2005], which gives us the
${\mathcal{O}}\left(\varepsilon^{-2(1+\delta)}\right)$ sampling complexity for
the LBL case for $0<\delta<1/2$.

######  Theorem A.4.

Suppose that $X,X_{1},X_{2},...$ are i.i.d. random variables., and set
$S_{n}=\sum_{k=1}^{n}X_{k},n\geq 1$. If $\mathbf{E}|X|^{p}<\infty$ and
$\mathbf{E}X=0$ when $1\leq p<2$, then

$\mathbf{E}\left|\frac{S_{n}}{n^{1/p}}\right|^{p}=\mathbf{E}\frac{|S_{n}|^{p}}{n}\stackrel{{\scriptstyle
n\to\infty}}{{\longrightarrow}}0.$

##  Appendix B Proof of Theorem 2.2

###### Proof.

Case 1: $d=D$
When $d=D$, Algorithm 1 samples one $y^{(D)}\sim\pi_{D}$ and outputs
$R_{D}(y^{(0:D-1)}):=g_{D}(y^{(0:D)})$. We first prove our output
$R_{D}(y^{(0:D-1)})$ has a finite expectation for almost surely every fixed
$y^{(0:D-1)}$, then its expectation equals $\gamma_{D}(y^{(0:D-1)})$ follows
directly from the algorithm design. To show the first point, notice that the
expectation of $|R_{D}(y^{(0:D-1)})|$ given $y^{(0:D-1)}$ equals the
conditional expectation $\mathbf{E}[\lvert g_{D}(y^{(0:D)})\rvert\mid
y^{(0:D-1)}]$. Since $\mathbf{E}[|g_{D}|]<\infty$ by assumption, we have
$\mathbf{E}[\lvert g_{D}(y^{(0:D)})\rvert\mid y^{(0:D-1)}]<\infty$ almost
surely. Therefore Algorithm 1 is unbiased when $d=D$ for almost surely every
input $y^{(0:D-1)}$. Furthermore, it has computational cost $1$, and the
output has finite $2^{D+1}$-st moment.

Case 2: $0\leq d\leq D-1$
Now that our base case is proven, we proceed via backwards induction. Suppose
(a), (b), (c) are satisfied for $d+1$ where $0\leq d\leq D-1$. Conditioning on
$y^{(0:d-1)}$, we sample $y^{(d)}\sim\pi_{d}$ and
$N_{d}\sim{\mathsf{Geo}}(r_{d})$. Algorithm 1 will call itself independently
for $2^{N_{d}}$ times, each with input {Depth index: $d+1$, Trajectory
History: $H=y^{(0:d)}$, Parameters: $r_{d+1},\cdots,r_{D-1}$}. This gives us
i.i.d. samples $R_{d+1}(y^{(0:d)})(1),...,R_{d+1}(y^{(0:d)})(2^{N_{d}})$ which
are used to compute the following:

$\displaystyle S_{2^{N_{d}}}$
$\displaystyle=R_{d+1}(y^{(0:d)})(1)+R_{d+1}(y^{(0:d)})(2)+\cdots+R_{d+1}(y^{(0:d)})(2^{N_{d}}),$
$\displaystyle S^{\mathsf{O}}_{2^{N_{d}-1}}$
$\displaystyle=R_{d+1}(y^{(0:d)})(1)+R_{d+1}(y^{(0:d)})(3)+\cdots+R_{d+1}(y^{(0:d)})(2^{N_{d}}-1),$
$\displaystyle S^{\mathsf{E}}_{2^{N_{d}-1}}$
$\displaystyle=R_{d+1}(y^{(0:d)})(2)+R_{d+1}(y^{(0:d)})(4)+\cdots+R_{d+1}(y^{(0:d)})(2^{N_{d}}).$

Then Algorithm 1 returns as output
$R_{d}(y^{(0:d)})=\Delta_{N_{d}}/p_{r_{d}}(N_{d})$, where $\Delta_{N_{d}}$ is
the antithetic quantity in Algorithm 1. By the inductive hypothesis on $d+1$,
we have for almost surely every $y^{(0:d)}$:

$\mathbf{E}_{\pi_{d+1:D}}[R_{d+1}(y^{(0:d)})\mid
y^{(0:d)}]=\gamma_{d+1}(y^{(0:d)})$

and

$\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]<\left(\prod_{i=d+1}^{D}\tilde{C}_{i}\right)\norm{g_{D}(y^{(0:D)})}_{\pi,2^{D+1}}^{2^{D+1}}.$

We will start with showing $R_{d}(y^{(0:d-1)})$ has a finite computational
cost and a finite $2^{d+1}$-th moment, and then show the unbiasedness.

Finite cost:
To show the computational cost, recall that implementing Algorithm 1 with
input depth $d$ requires $2^{N_{d}}$ calls of Algorithm 1 with input depth
$d+1$. Since $N_{d}\sim{\mathsf{Geo}}(r_{d})$ with $r_{d}>0.5$, calling
Algorithm 1 with input depth $d$ has an expected cost:

$\displaystyle\frac{r_{d}}{2r_{d}-1}\times\text{the expected cost of Algorithm
\ref{alg:recursive-rMLMC} with input depth}~{}d+1,$

where $\frac{r_{d}}{2r_{d}-1}=\mathbf{E}[2^{N_{d}}]<\infty$. By our inductive
hypothesis, the second term in the above product is finite, therefore the
expected cost of Algorithm 1 with input depth $d$ is also finite.

Finite $2^{d+1}$-th moment:

Next we show $R_{d+1}(y^{(0:d)})$ has a finite $2^{d+1}$-th moment. For every
fixed positive integer $n$, doing a Taylor expansions for $g_{d}$ at
$(y^{(0:d)},\gamma_{d+1})$ with respect to the last component gives us:

$\displaystyle g_{d}\left(y^{(0:d)},\frac{S_{2^{n}}}{2^{n}}\right)$
$\displaystyle=g_{d}(y^{(0:d)},\gamma_{d+1})+\partial_{d+2}g_{d}(y^{(0:d)},\gamma_{d+1})\left(\frac{S_{2^{n}}}{2^{n}}-\gamma_{d+1}\right)$
$\displaystyle+\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi(n))\left(\frac{S_{2^{n}}}{2^{n}}-\gamma_{d+1}\right)^{2}$
$\displaystyle
g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}\right)$
$\displaystyle=g_{d}(y^{(0:d)},\gamma_{d+1})+\partial_{d+2}g_{d}(y^{(0:d)},\gamma_{d+1})\left(\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)$
$\displaystyle\quad+\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi^{\mathsf{O}}(n-1))\left(\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)^{2}$
$\displaystyle
g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}\right)$
$\displaystyle=g_{d}(y^{(0:d)},\gamma_{d+1})+\partial_{d+2}g_{d}(y^{(0:d)},\gamma_{d+1})\left(\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)$
$\displaystyle\quad+\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi^{\mathsf{E}}(n-1))\left(\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)^{2},$

with $\xi(n)$ between $\gamma_{d+1}$ and $S_{2^{n}}/2^{n}$,
$\xi^{\mathsf{O}}(n-1)$ between $\gamma_{d+1}$ and
$S^{\mathsf{O}}_{2^{n-1}}/2^{n-1}$, $\xi^{\mathsf{E}}(n-1)$ between
$\gamma_{d+1}$ and $S_{2^{n-1}}/2^{n-1}$.

Thus, we have:

$\begin{split}\Delta_{n}&={}g_{d}\left(y^{(0:d)},\frac{S_{2^{n}}}{2^{n}}\right)-\frac{1}{2}\left[g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{O}}_{2^{{n}-1}}}{2^{{n}-1}}\right)+g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{E}}_{2^{{n}-1}}}{2^{{n}-1}}\right)\right]\\\
&=\partial_{d+2}g_{d}(y^{(0:d)},\gamma_{d+1})\left(\frac{S_{2^{n}}}{2^{n}}-\gamma_{d+1}\right)+\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi(n))\left(\frac{S_{2^{n}}}{2^{n}}-\gamma_{d+1}\right)^{2}\\\
&\quad-\frac{1}{2}\Bigg{[}\partial_{d+2}g_{d}(y^{(0:d)},\gamma_{d+1})\left(\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)+\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi^{\mathsf{O}}(n-1))\left(\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)^{2}\\\
&\quad\quad\quad+\partial_{d+2}g_{d}(y^{(0:d)},\gamma_{d+1})\left(\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)+\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi^{\mathsf{E}}(n-1))\left(\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)^{2}\Bigg{]}\\\
&=\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi(n))\left(\frac{S_{2^{n}}}{2^{n}}-\gamma_{d+1}\right)^{2}-\frac{1}{2}\Bigg{[}\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi^{\mathsf{O}}(n-1))\left(\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)^{2}\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{1}{2}\partial_{d+2}^{2}g_{d}(y^{(0:d)},\xi^{\mathsf{E}}(n-1))\left(\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)^{2}\Bigg{]}.\end{split}$

By the LBS assumption which assumes
$\lvert\partial_{d+2}^{2}g_{d}(y^{(0:d)},z))\rvert<K_{d}$ for every
$(y^{(0:d)},z)$, and our inductive hypothesis:
$\gamma_{d+1}=\mathbf{E}[R_{d+1}(y^{(0:d)})\mid y^{(0:d)}]$ which allows us to
use Lemma A.3 with $Z_{1}=y^{(0:d)}$, $Z_{2}=R_{d+1}(y^{(0:d)})$, we have:

$\displaystyle\norm{\Delta_{n}}_{\pi,2^{d+1}}$ $\displaystyle\leq
K_{d}\norm{\left(\frac{S_{2^{n}}}{2^{n}}-\gamma_{d+1}\right)^{2}}_{\pi,2^{d+1}}+K_{d}\norm{\left(\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}\right)^{2}}_{\pi,2^{d+1}}$
$\displaystyle\leq
K_{d}B^{\prime}_{2^{d+2}}\left(\frac{\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]}{2^{(2^{d+1}n)}}\right)^{1/2^{d+1}}+K_{d}B^{\prime}_{2^{d+2}}\left(\frac{\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]}{2^{(2^{d+1}(n-1))}}\right)^{1/2^{d+1}}$
$\displaystyle=K_{d}B^{\prime}_{2^{d+2}}\left(\frac{\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]}{2^{(2^{d+1}n)}}\right)^{1/2^{d+1}}+2K_{d}B^{\prime}_{2^{d+2}}\left(\frac{\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]}{2^{(2^{d+1}n)}}\right)^{1/2^{d+1}}$
$\displaystyle=3K_{d}B^{\prime}_{2^{d+2}}\left(\frac{\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]}{2^{(2^{d+1}n)}}\right)^{1/2^{d+1}}.$

Therefore, in total:

$\mathbf{E}_{\pi}\left[|\Delta_{n}|^{2^{d+1}}\right]\leq
D_{d}\frac{\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]}{2^{(2^{d+1}n)}},$

with $D_{d}=(3K_{d}B^{\prime}_{2^{d+2}})^{(2^{d+1})}$.

The result claimed in (c) is obtained as follows. We have for
$(1-r_{d})=2^{-k_{d}}$ for some
$k_{d}\in\left(1,\frac{2^{d+1}}{2^{d+1}-1}\right)$:

$\displaystyle\mathbf{E}_{\pi}\left[|R_{d}(y^{(0:d-1)})|^{2^{d+1}}\right]$
$\displaystyle=\sum_{n=0}^{\infty}\frac{\mathbf{E}_{\pi}\left[|\Delta_{n}|^{2^{d+1}}\right]}{p_{r_{d}}(n)^{{2^{d+1}}-1}}$
$\displaystyle\leq\frac{D_{d}\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2^{d+2}}\right]}{r_{d}^{2^{d+1}-1}}\sum_{n=0}^{\infty}\frac{1}{2^{2^{d+1}n}(1-r_{d})^{({2^{d+1}-1})n}}$
$\displaystyle\leq
C_{d}\left(\prod_{i=d+1}^{D}\tilde{C}_{i}\right)\norm{g_{D}^{2}(y^{(0:D)})}_{\pi,2^{D}}^{2^{D}}\sum_{n=0}^{\infty}\left(\frac{1}{2^{{2^{d+1}}-k_{d}(2^{d+1}-1)}}\right)^{n}$
$\displaystyle=C_{d}\left(\prod_{i=d+1}^{D}\tilde{C}_{i}\right)\norm{g_{D}^{2}(y^{(0:D)})}_{\pi,2^{D}}^{2^{D}}\left(\frac{2^{({2^{d+1}}-k_{d}(2^{d+1}-1)})}{2^{({2^{d+1}}-k_{d}(2^{d+1}-1))}-1}\right)$
$\displaystyle=\left(\prod_{i=d}^{D}\tilde{C}_{i}\right)\norm{g_{D}^{2}(y^{(0:D)})}_{\pi,2^{D}}^{2^{D}},$

here the second inequality follows from our inductive hypothesis. The choice
$k_{d}\in\left(1,\frac{2^{d+1}}{2^{d+1}-1}\right)$ is crucial. It ensures
$2^{(2^{d+1})}(1-r_{d})^{(2^{d+1}-1)}>1$, and in turn ensures the infinite
summation of the above geometric series is finite.

Unbiasedness:
Now we show the unbiasedness of $R_{d}(y^{(0:d-1)})$. Firstly, since we have
just shown $R_{d}(y^{(0:d-1)})$ has a finite $2^{(}d+1)$-th moment under
$\pi$, it directly implies $\lvert R_{d}(y^{(0:d-1)})\rvert$ has a finite
first moment, which further implies

$\mathbf{E}[\lvert R_{d}(y^{(0:d-1)})\rvert\mid y^{(0:d-1)}]$

is finite for $\pi$-almost surely $y^{(0:d-1)}$.

We fix $y^{(0:d-1)}$ from now on, and we will write
$\mathbf{E}_{\pi_{d:D}}[\cdot]$ as a shorthand notation for
$\mathbf{E}[\cdot\mid y^{(0:d-1)}]$ . Recall that we have output
$R_{d}=\Delta_{N_{d}}/p_{r_{d}}(N_{d})$, with

$\Delta_{N_{d}}=g_{d}\left(y^{(0:d)},\frac{S_{2^{N_{d}}}}{2^{N_{d}}}\right)-\frac{1}{2}\left[g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{O}}_{2^{{N_{d}}-1}}}{2^{{N_{d}}-1}}\right)+g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{E}}_{2^{{N_{d}}-1}}}{2^{{N_{d}}-1}}\right)\right].$

Then,

$\displaystyle\mathbf{E}_{\pi_{d:D}}[R_{d}(y^{(0:d-1)})]$
$\displaystyle=\mathbf{E}_{\pi_{d:D}}\left[\mathbf{E}\left[\frac{\Delta_{N_{d}}}{p_{r_{d}}(N_{d})}\mid
N_{d}\right]\right]$
$\displaystyle=\mathbf{E}_{\pi_{d:D}}\left[\sum_{n=0}^{\infty}\frac{\Delta_{n}}{p_{r_{d}}(n)}p_{r_{d}}(n)\right]$
$\displaystyle=\mathbf{E}_{\pi_{d:D}}\left[\sum_{n=0}^{\infty}\Delta_{n}\right]$
$\displaystyle\stackrel{{\scriptstyle(\star\star\star)}}{{=}}\sum_{n=0}^{\infty}\mathbf{E}_{\pi_{d:D}}[\Delta_{n}]$
$\displaystyle=\sum_{n=0}^{\infty}\mathbf{E}_{\pi_{d:D}}\left\\{g_{d}\left(y^{(0:d)},\frac{S_{2^{n}}}{2^{n}}\right)-\frac{1}{2}\left[g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{O}}_{2^{{n}-1}}}{2^{{n}-1}}\right)+g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{E}}_{2^{{n}-1}}}{2^{{n}-1}}\right)\right]\right\\}$
$\displaystyle=\mathbf{E}_{\pi_{d:D}}\left[g_{d}\left(y^{(0:d)},\lim_{n\to\infty}\frac{S_{2^{n}}}{2^{n}}\right)\right]$
$\displaystyle=\mathbf{E}_{\pi_{d:D}}\left[g_{d}\left(y^{(0:d)},\gamma_{d+1}(y^{(0:d)})\right)\right]$
$\displaystyle=\gamma_{d}(y^{(0:d-1)}).$

All the above calculations are straightforward except for $(\star\star\star)$,
which swaps the order of expectation and summation. Therefore we complete this
proof of unbiasedness by justifying the swap in $(\star\star\star)$. To
justify the swap, it suffices to show
$\sum_{n}\mathbf{E}[|\Delta_{n}|]<\infty$. Notice that $R_{d}(y^{(0:d)})$ can
be equivalently written as
$\sum_{n=1}^{\infty}\Delta_{n}\mathbf{I}(N_{d}=n)/p_{r_{d}}(n)$ where $N_{d}$
independent with $\\{\Delta_{i}\\}$. Calculating
$\mathbf{E}_{\pi_{d:D}}[\lvert R_{d}(y^{(0:d)})\rvert]$ yields:

$\displaystyle\mathbf{E}_{\pi_{d:D}}[\lvert R_{d}(y^{(0:d)})\rvert]$
$\displaystyle=\mathbf{E}_{\pi_{d:D}}\left[\left\lvert\sum_{n=1}^{\infty}\frac{\Delta_{n}\mathbf{I}(N_{d}=n)}{p_{r_{d}}(n)}\right\rvert\right]$
$\displaystyle=\mathbf{E}_{\pi_{d:D}}\left[\sum_{n=1}^{\infty}\left\lvert\frac{\Delta_{n}\mathbf{I}(N_{d}=n)}{p_{r_{d}}(n)}\right\rvert\right]\quad{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{only
one term in the summation is non-zero}}$
$\displaystyle=\sum_{n=1}^{\infty}\mathbf{E}_{\pi_{d:D}}\left[\left\lvert\frac{\Delta_{n}\mathbf{I}(N_{d}=n)}{p_{r_{d}}(n)}\right\rvert\right]\quad{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{every
term is non-negative}}$
$\displaystyle=\sum_{n=1}^{\infty}\mathbf{E}_{\pi_{d:D}}\left[\left\lvert\frac{\Delta_{n}}{p_{r_{d}}(n)}\right\rvert\right]\mathbf{E}\left[\mathbf{I}(N_{d}=n)\right]\quad{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{independence
between $N$ and $\\{\Delta_{i}\\}$}}$
$\displaystyle=\sum_{n=1}^{\infty}\mathbf{E}_{\pi_{d:D}}\left[\left\lvert\Delta_{n}\right\rvert\right].$

Since we already know $\mathbf{E}_{\pi_{d:D}}[\lvert
R_{d}(y^{(0:d-1)})\rvert]<\infty$, this justifies our swap
$(\star\star\star)$.

∎

##  Appendix C Proof of Theorem 2.4

The proof strategy of Theorem 2.4 is very similar to Theorem 2.2. We start
with a backward induction.

###### Proof.

Case 1: $d=D$

When $d=D$, Algorithm 1 samples one $y^{(D)}\sim\pi_{D}$ and outputs
$R_{D}(y^{(0:D-1)}):=g_{D}(y^{(0:D)})$. Again, we first prove our output
$R_{D}(y^{(0:D-1)})$ has a finite expectation for almost surely every fixed
$y^{(0:D-1)}$, then its expectation equals $\gamma_{D}(y^{(0:D-1)})$ follows
directly from the algorithm design. To show the first point, notice that the
expectation of $|R_{D}(y^{(0:D-1)})|$ equals the conditional expectation
$\mathbf{E}[\lvert g_{D}(y^{(0:D)})\rvert\mid y^{(0:D-1)}]$. Since
$\mathbf{E}[|g_{D}|]<\infty$ by assumption, we have $\mathbf{E}[\lvert
g_{D}(y^{(0:D)})\rvert\mid y^{(0:D-1)}]<\infty$ almost surely. Therefore
Algorithm 1 is unbiased when $d=D$ for almost surely every input
$y^{(0:D-1)}$. Furthermore, it has computational cost $1$, and the output has
finite $\left(2-\frac{\delta}{2^{D}}\right)$-th moment.

Case 2: $0\leq d\leq D-1$

Now that our base case is proven, we proceed via backwards induction. Let
$\delta_{d}:=\delta/2^{d}$ for every $d\in\\{0,1,\ldots,D\\}$. Suppose
unbiasedness, finite $(2-\delta_{d+1})$-th moment, and finite expected
computational cost are all satisfied for $d+1$ where $0\leq d\leq D-1$. Then
Algorithm 1 will call itself independently for $2^{N_{d}}$ times, each with
input {Depth index: $d+1$, Trajectory History: $H=y^{(0:d)}$, Parameters:
$r_{d+1},\cdots,r_{D-1}$}. This gives us i.i.d. samples
$R_{d+1}(y^{(0:d)})(1),...,R_{d+1}(y^{(0:d)})(2^{N_{d}})$ which are used to
compute the following:

$\displaystyle S_{2^{N_{d}}}$
$\displaystyle=R_{d+1}(y^{(0:d)})(1)+R_{d+1}(y^{(0:d)})(2)+\cdots+R_{d+1}(y^{(0:d)})(2^{N_{d}}),$
$\displaystyle S^{\mathsf{O}}_{2^{N_{d}-1}}$
$\displaystyle=R_{d+1}(y^{(0:d)})(1)+R_{d+1}(y^{(0:d)})(3)+\cdots+R_{d+1}(y^{(0:d)})(2^{N_{d}}-1),$
$\displaystyle S^{\mathsf{E}}_{2^{N_{d}-1}}$
$\displaystyle=R_{d+1}(y^{(0:d)})(2)+R_{d+1}(y^{(0:d)})(4)+\cdots+R_{d+1}(y^{(0:d)})(2^{N_{d}}).$

Then Algorithm 1 returns as output
$R_{d}(y^{(0:d)})=\Delta_{N_{d}}/p_{r_{d}}(N_{d})$, where $\Delta_{N_{d}}$ is
defined in Algorithm 1. By the inductive hypothesis on $d+1$, we have for
almost surely every $y^{(0:d)}$:

$\mathbf{E}[R_{d+1}(y^{(0:d)})\mid y^{(0:d)}]=\gamma_{d+1}(y^{(0:d)})$

and

$\mathbf{E}_{\pi}\left[|R_{d+1}(y^{(0:d)})|^{2-\delta_{d+1}}\right]<\left(\prod_{i=d}^{D}\tilde{C}_{i}\right)\norm{g_{D}(y^{(0:D)})}_{\pi,2}^{2-\delta_{d+1}}.$

We will start with showing $R_{d}(y^{(0:d-1)})$ has a finite computational
cost and a finite $(2-\delta_{d})$-th moment, and then show the unbiasedness.
Finite cost:
To show the computational cost, recall that implementing Algorithm 1 with
input depth $d$ requires $2^{N_{d}}$ calls of Algorithm 1 with input depth
$d+1$. It suffices to check $r_{d}>0.5$, which reduces to check the the upper
bound for $k_{d}$ (defined in Theorem 2.4) staisfies

$\left(\frac{2^{d+2}-3\delta}{2^{d+3}-3\delta}\right)\left(\frac{2^{d+1}-\delta}{2^{d}-\delta}\right)>1.$

Let $t\coloneqq 2^{d}$ and the above product becomes:

$\frac{4t-3\delta}{8t-3\delta}\frac{2t-\delta}{t-\delta}=\frac{8t^{2}+3\delta^{2}-10\delta}{8t^{2}+3\delta^{2}-11\delta}>1.$

Since $N_{d}\sim{\mathsf{Geo}}(r_{d})$ with $r_{d}>0.5$, calling Algorithm 1
with input depth $d$ has an expected cost:

$\displaystyle\frac{r_{d}}{2r_{d}-1}\times\text{the expected cost of Algorithm
\ref{alg:recursive-rMLMC} with input depth}~{}d+1,$

where $\frac{r_{d}}{2r_{d}-1}=\mathbf{E}[2^{N_{d}}]<\infty$. By our inductive
hypothesis, the second term in the above product is finite, therefore the
expected cost of Algorithm 1 with input depth $d$ is also finite.

Finite $(2-\delta_{d})$-th moment:
Next we show $R_{d}$ has a finite $(2-\delta_{d})$-th moment. By the uniform
$L_{d}$-Lipschitz property of $g_{d}$:

$\displaystyle|\Delta_{n}|$
$\displaystyle\leq\frac{1}{2}\left|g_{d}\left(y^{(0:d)},\frac{S_{2^{n}}}{2^{n}}\right)-g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}\right)\right|+\frac{1}{2}\left|g_{d}\left(y^{(0:d)},\frac{S_{2^{n}}}{2^{n}}\right)-g_{d}\left(y^{(0:d)},\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}\right)\right|$
$\displaystyle\leq\frac{L_{d}}{2}\left|\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}\right|.$

Therefore, for any fixed $1\leq p<2$, applying triangle inequality under the
norm $\lVert\cdot\rVert_{\pi,p}$, and applying Lemma A.3 with
$Z_{1}=y^{(0:d)}$, $Z_{2}=R_{d+1}(y^{(0:d)})$, we have:

$\displaystyle\norm{\Delta_{n}}_{\pi,p}$
$\displaystyle\leq\frac{L_{d}}{2}\norm{\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}}_{\pi,p}$
$\displaystyle\leq\frac{L_{d}}{2}\norm{\frac{S^{\mathsf{O}}_{2^{n-1}}}{2^{n-1}}-\gamma_{d+1}}_{\pi,p}+\frac{L_{d}}{2}\norm{\gamma_{d+1}-\frac{S^{\mathsf{E}}_{2^{n-1}}}{2^{n-1}}}_{\pi,p}$
$\displaystyle\leq
L_{d}\left(\frac{B_{p}\mathbf{E}_{\pi}[|R_{d+1}(y^{(0:d)})|^{p}]}{2^{(n-1)(p-1)}}\right)^{1/p},$

exponentiate both sides by $p$ yields,

$\mathbf{E}_{\pi}[|\Delta_{n}|^{p}]\leq\frac{L_{d}^{p}B_{p}\mathbf{E}_{\pi}[|R_{d+1}(y^{(0:d)})|^{p}]}{2^{(n-1)(p-1)}}\leq\frac{C(d,p)\mathbf{E}_{\pi}[|R_{d+1}(y^{(0:d)})|^{p}]}{2^{(p-1)n}},$

where $C(d,p)=L_{d}^{p}B_{p}2^{1-p}$.

Recall that $\delta_{d}=\delta/2^{d}$, let us choose
$q_{d}=2-(\delta_{d}+\delta_{d+1})/2$. Since

$\left(1,\left(\frac{q_{d}-1}{q_{d}}\right)\left(\frac{2-\delta_{d}}{1-\delta_{d}}\right)\right)=\left(1,\left(\frac{2^{d+2}-3\delta}{2^{d+3}-3\delta}\right)\left(\frac{2^{d+1}-\delta}{2^{d}-\delta}\right)\right),$

by definition of $k_{d}$ in the Theorem statement we have

$k_{d}<\left(\frac{q_{d}-1}{q_{d}}\right)\left(\frac{2-\delta_{d}}{1-\delta_{d}}\right).$

Now we estimate the $(2-\delta_{d})$-th moment of $R_{d}$. An important trick
in the calculation below is that we are not going to use the above estimate of
$\mathbf{E}_{\pi}[\lvert\Delta_{n}\rvert^{p}]$ directly on $p=2-\delta_{d}$.
Instead, we will first use Hölder’s inequality, and then bound
$\mathbf{E}_{\pi}[|\Delta_{n}|^{q_{d}}]^{(2-\delta_{d})/q_{d}}$ via the above
estimate. It turns out the first way gives us an order of
$2^{-n(1-\delta_{d})}$, while the latter is of order
$2^{-n(q_{d}-1)(2-\delta_{d})/q_{d}}$. Since the function
$(x-1)(2-\delta_{d})/x$ is increasing with $x$ when $x>1$, and equals
$1-\delta_{d}$ when $x=2-\delta_{d}$, we gain an extra factor
$2^{-\Omega(1)n}$ by choosing $q_{d}>2-\delta_{d}$ and use Hölder’s
inequality, which is important for establishing our main result. The detailed
calculation is below:

$\displaystyle\mathbf{E}_{\pi}$
$\displaystyle[|R_{d}(y^{(0:d-1)})|^{2-\delta_{d}}]\leq\sum_{n=0}^{\infty}\frac{\mathbf{E}_{\pi}[|\Delta_{n}|^{2-\delta_{d}}]}{p_{r_{d}}(n)^{1-\delta_{d}}}$
$\displaystyle\leq\sum_{n=0}^{\infty}\frac{\left(\mathbf{E}_{\pi}\left[|\Delta_{n}|^{2-\delta_{d}\cdot\frac{q_{d}}{2-\delta_{d}}}\right]\right)^{\frac{2-\delta_{d}}{q_{d}}}}{p_{r_{d}}(n)^{1-\delta_{d}}}\quad\quad\quad\quad\quad\quad{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{H\"{o}lder's
inequality}}$
$\displaystyle\leq\frac{1}{r_{d}^{1-\delta_{d}}}\sum_{n=0}^{\infty}\left(\frac{C(d,q_{d})\mathbf{E}_{\pi}[|R_{d+1}(y^{(0:d)})|^{q_{d}}]}{2^{(q_{d}-1)n}}\right)^{\frac{2-\delta_{d}}{q_{d}}}\frac{1}{(1-r_{d})^{(1-\delta_{d})n}}\quad{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{estimate
of $\mathbf{E}_{\pi}[|\Delta_{n}|^{p}]$}}$
$\displaystyle=C^{\prime}(d)\norm{R_{d+1}(y^{(0:d)})}^{2-\delta_{d}}_{\pi,q_{d}}\sum_{n=0}^{\infty}\left(\frac{1}{2^{\frac{(q_{d}-1)}{q_{d}}(2-\delta_{d})-k_{d}(1-\delta_{d})}}\right)^{n}\quad\qquad~{}~{}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{here}~{}C^{\prime}(d)=\frac{C(d,q_{d})^{(2-\delta_{d})/q_{d}}}{r_{d}^{1-\delta_{d}}}}$
$\displaystyle\leq
C^{\prime}(d)\norm{R_{d+1}(y^{(0:d)})}^{2-\delta_{d}}_{\pi,2-\delta_{d+1}}\sum_{n=0}^{\infty}\left(\frac{1}{2^{\frac{(q_{d}-1)}{q_{d}}(2-\delta_{d})-k_{d}(1-\delta_{d})}}\right)^{n}\quad\quad~{}~{}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{since}~{}q_{d}<2-\delta_{d+1}}$
$\displaystyle\leq
C^{\prime}(d)\left(\prod_{i=d+1}^{D}\tilde{C}_{i}\right)\norm{g_{D}(y^{(0:D)})}_{\pi,2}^{2-\delta_{d}}\left(\frac{2^{\frac{(q_{d}-1)}{q_{d}}(2-\delta_{d})-k_{d}(1-\delta_{d})}}{2^{\frac{(q_{d}-1)}{q_{d}}(2-\delta_{d})-k_{d}(1-\delta_{d})}-1}\right)\quad{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\text{inductive
hypothesis}}$
$\displaystyle=\left(\prod_{i=d}^{D}\tilde{C}_{i}\right)\norm{g_{D}(y^{(0:D)})}_{\pi,2}^{2-\delta_{d}},$

and note the RHS is still finite given the assumption of our theorem on
$g_{D}$. Again, as we can see in the proof, the choice of $k_{d}$ and $q_{d}$
is crucial for our calculation. It ensures
$\frac{(q_{d}-1)}{q_{d}}(2-\delta_{d})-k_{d}(1-\delta_{d})>0$, and in turn
ensures the above summation of the geometric series converges.

Unbiasedness:

The proof of unbiasedness of our estimator in this case is identical to the
LBS case, however we still require a justification of the existence of a
finite conditional expectation of $R_{d}(y^{(0:d-1)})$. By what we have just
proven,

$\mathbf{E}_{\pi}[|R_{d}(y^{(0:d-1)})|]\leq\left(\mathbf{E}_{\pi}\left[|R_{d}(y^{(0:d-1)})|^{2-\frac{\delta}{2^{d}}}\right]\right)^{1/(2-\frac{\delta}{2^{d}})}<\infty.$

Given $\mathbf{E}_{\pi}[|R_{d}(y^{(0:d-1)})|]<\infty$, we immediately have
$\mathbf{E}_{\pi_{d:D}}[R_{d}(y^{(0:d-1)})]$ exists for almost surely every
$y^{(0:d-1)}$.

∎

##  Appendix D Construction of the NMC estimator

The construction of the NMC estimator is described in Rainforth et al. [2018].
For concreteness, we explain the construction details here for the $D=2$ case,
which we use in Section 3.

Fix positive integers $N_{0},N_{1},N_{2}$, users first simulate $N_{0}$ i.i.d.
$\\{y^{(0)}_{i}\\}_{i=1}^{N_{0}}\sim\pi_{0}$. For each fixed $y^{(0)}_{i}$,
users sample $N_{1}$ i.i.d. $\\{y^{(1)}_{i,j}\\}_{j=1}^{N_{1}}$ from
$\pi_{1}(\cdot\mid y^{(0)}_{i})$. Then for each fixed trajectory
$(y^{(0)}_{i},y^{(1)}_{i,j})$, users sample $N_{2}$ i.i.d.
$\\{y^{(2)}_{i,j,k}\\}_{k=1}^{N_{2}}$ from $\pi_{2}(\cdot\mid
y^{(0)}_{i},y^{(1)}_{i,j})$. After getting all these samples, we use the
standard estimator:

$\hat{\gamma}_{2}(y^{(0)}_{i},y^{(1)}_{i,j})\coloneqq\frac{1}{N_{2}}\sum_{k=1}^{N_{2}}g_{2}\left(y^{(0)}_{i},y^{(1)}_{i,j},y^{(2)}_{i,j,k}\right)$

to estimate $\gamma_{2}(y^{(0)}_{i},y^{(1)}_{i,j})$. Then using the plug-in
estimator

$\hat{\gamma}_{1}(y^{(0)}_{i})\coloneqq\frac{1}{N_{1}}\sum_{j=1}^{N_{1}}g_{1}\left(y^{(0)}_{i},y^{(1)}_{i,j},\hat{\gamma}_{2}(y^{(0)}_{i},y^{(1)}_{i,j})\right)$

to estimate $\gamma_{1}(y^{(0)}_{i})$. Finally, plugging in these $N_{0}$
estimators $\\{\hat{\gamma}_{1}(y^{(0)}_{i})\\}_{i=1}^{N_{0}}$ to form the NMC
estimator

$\hat{\gamma}_{0}\coloneqq\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}g_{0}\left(y^{(0)}_{i},\hat{\gamma}_{1}(y^{(0)}_{i})\right)$

for $\gamma_{0}$.

It is proven in Rainforth et al. [2018] that when all $N_{0},N_{1},N_{2}$ go
to $\infty$, $\hat{\gamma}_{0}$ converges to $\gamma_{0}$. It remains crucial
to allocate $N_{0},N_{1},N_{2}$ to maximize the convergence rate with respect
to the total sample size $n=N_{0}N_{1}N_{2}$. The choice
$N_{0}=N_{1}^{2}=N_{2}^{2}$ is suggested in Rainforth et al. [2018], which has
a ${\mathcal{O}}(N^{-1/4})$ convergence rate for the rMSE, or a
${\mathcal{O}}(N^{-1/2})$ rate for the MSE.

##  Appendix E Additional statistics of Section 3

Figure 1 in Section 3 compares the errors between $\mathsf{READ}$, NMC1, and
NMC2 in terms of the cost of the total sample size. Here we also compare their
estimation errors in terms of the wall-clock time. For $\mathsf{READ}$, we
call Algorithm 1 repeatedly for $10^{5}$ times. For NMC1, we choose
$N_{0}=N_{1}=N_{2}=400$. For NMC2, we choose
$N_{0}=N_{1}^{2}=N_{2}^{2}=10^{4}$. Their estimation errors and the
corresponding wall-clock time costs are summarized in Table 1. It is clear
that $\mathsf{READ}$ is both faster (in wall-clock time) and more accurate. We
also calculate the time-normalized squared error, defined as the product
between the time cost and the squared error Glynn and Whitt [1992]. From the
normalized squared error, $\mathsf{READ}$ is more than $130$ times more
efficient than NMC2, and more than $407$ times more efficient than NMC1.

Method | Setting |  | Total Sample
---
Cost
Time/Seconds | Squared Error |  | Time-normalized
---
Squared Error
$\mathsf{READ}$ |  $10^{5}$ repetitions | $4.625\times 10^{5}$ | $\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{11.86}$ | $\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{3.186\times 10^{-6}}$ | $\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{3.78\times 10^{-5}}$
NMC1 | $N_{0}=N_{1}=N_{2}=400$ | $6.4\times 10^{7}$ | $43.79$ | $3.51\times 10^{-4}$ | $1.54\times 10^{-2}$
NMC2 | $N_{0}=N_{1}^{2}=N_{2}^{2}=10^{4}$ | $10^{8}$ | $72.6$ | $6.8\times 10^{-5}$ | $4.93\times 10^{-3}$

Table 1: Cost comparison between different methods

##  Appendix F Extra experiments

We consider an extra experiment with unknown ground truth. Let
$\sigma(x)\coloneqq e^{x}/(1+e^{x})$ be the sigmoid function. Suppose the
process $(y^{(0)},y^{(1)},y^{(2)})$ satisfies
$y^{(0)}\sim\mathbf{N}(0,1),y^{(1)}\sim\mathbf{N}(y^{(0)},1),y^{(2)}\sim\mathbf{N}(y^{(1)},1)$.
Define
$g_{0}(y^{(0)},z)\coloneqq\sigma(y^{(0)}+z),g_{1}(y^{(0:1)},z)\coloneqq\sigma(y^{(1)}+z)$,
and $g_{2}(y^{(0:2)})\coloneqq\sigma(y^{(2)})$. The target quantity
$\gamma_{0}$ defined (1) is again a nested expectation with $D=2$. Although we
can not analytically calculate out $\gamma_{0}$, we still implement our
$\mathsf{READ}$ estimator with the NMC estimators in Rainforth et al. [2018]
and compare their performance. The parameters of $\mathsf{READ}$ are the same
as Section 3, the allocation of $N_{0},N_{1},N_{2}$ of the NMC estimators also
follows the same way as Section 3.

The scatter plot of the estimation results is shown in Figure 3. Although no
ground truth is available, the trend for the estimation is clear. All three
methods eventually get close to $0.612$, represented by the dotted black line
in Figure 3. It is also clear from the plot that $\mathsf{READ}$ always stays
very close to the black line. In contrast, both NMC1 and NMC2 are
significantly more fluctuated than $\mathsf{READ}$, where NMC1 appears to be
the most unstable estimator. This again matches with the theoretical
predictions in our paper and Rainforth [2018] that $\mathsf{READ}$ converges
the fastest while NMC1 converges the slowest.

Figure 3: Scatterplot of the estimation of $\gamma_{0}$ as a function of
$\log_{10}(\text{Total Sample Cost})$. Blue, red, green points correspond to
$\mathsf{READ}$, NMC$1$, NMC$2$ estimators respectively.

Figure 4: (a): Heatmap of the (empirical) standard deviation of
$\mathsf{READ}$. (b): Heatmap of the work-normalized standard deviation of
$\mathsf{READ}$. Here $r_{0}\in(0.6,0.74)$, $r_{1}\in(0.55,0.6)$. Each
standard deviation is estimated based on $10^{6}$ repetitions of Algorithm 1.

Next we let the parameters $(r_{0},r_{1})$ in Algorithm 1 vary and investigate
the proper choice of the parameters in this experiment. Theorem 2.2 shows any
$(r_{0},r_{1})\in(0.5,0.75)\times(0.5,1-2^{-4/3})$ guarantees $\mathsf{READ}$
has finite variance and finite cost. Therefore we choose $r_{0}$ on the
lattice $\\{0.6,0.614,\ldots,0.74\\}$ and $r_{1}$ on the lattice
$\\{0.55,0.555,\ldots,0.6\\}$. For each pair of $(r_{0},r_{1})$, we repeat
Algorithm 1 for $10^{6}$ times, record the results and calculate their
empirical standard deviation. The heatmap is shown in the left plot of Figure
4. The pattern suggests the standard deviation depends crucially on the choice
of $r_{0}$, but less on $r_{1}$. The standard deviation decreases when $r_{0}$
increases.

Since the expected sample cost of Algorithm 1 equals
$\left(r_{1}/(2r_{1}-1)\right)\left(r_{2}/(2r_{2}-1)\right)$. We also plot the
‘work-normalized standard deviation’, which is defined as
$\sqrt{\text{Expected Sample Cost}}\times\text{Standard deviation}$ in Glynn
and Whitt [1992] to measure the efficiency of difference choices of
$(r_{0},r_{1})$. The heatmap is shown in the right subplot of Figure 4. Our
result suggests users should choose larger values of $(r_{0},r_{1})$ to
maximize the efficiency, at least in this example.

Finally we test our results when $r_{0},r_{1}$ are both beyond the range given
by Theorem 2.2. Algorithm 1 can still be implemented, though there is no
guarantees on the finite variance. Nevertheless, we choose
$r_{0}\in\\{0.8,0.81,\ldots,0.9\\}$ and $r_{1}\in\\{0.7,0.71,\ldots,0.8\\}$
and report the heatmaps of the standard deviations/work-normalized standard
deviations in Figure s5. The estimates become significantly less stable, as
some pairs of $(r_{0},r_{1})$ have much larger standard deviation than their
neighborhoods. This suggests the actual standard deviation maybe already
infinity (though the empirical standard deviation will always be finite), and
therefore our result is less reliable. In conclusion, although larger values
of $(r_{0},r_{1})$ can reduce the average cost of each implementation, users
should not choose them too large as it may sacrifice the finite variance.
Users can choose the parameters closer to the upper end of the ranges in
Theorem 2.2, but not exceed these ranges.

Figure 5: (a): Heatmap of the (empirical) standard deviation of
$\mathsf{READ}$. (b): Heatmap of the work-normalized standard deviation of
$\mathsf{READ}$. Here $r_{0}\in(0.8,0.9)$, $r_{1}\in(0.7,0.8)$. Each standard
deviation is estimated based on $10^{6}$ repetitions of Algorithm 1.