# Towards a turnkey approach for unbiased Monte Carlo estimation of smooth functions of expectations

Nicolas Chopin, Francesca Crucinio, Sumeetpal Singh

► To cite this version:

Nicolas Chopin, Francesca Crucinio, Sumeetpal Singh. Towards a turnkey approach for unbiased Monte Carlo estimation of smooth functions of expectations. Biometrika, 2025, 112 (3), 10.1093/biomet/asaf030. hal-05365655

HAL Id: hal-05365655

https://hal.science/hal-05365655v1

Submitted on 14 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# TOWARDS A TURNKEY APPROACH TO UNBIASED MONTE CARLO ESTIMATION OF SMOOTH FUNCTIONS OF EXPECTATIONS

NICOLAS CHOPIN, FRANCESCA R. CRUCINIO, AND SUMEETPAL S. SINGH

**ABSTRACT.** Given a smooth function $f$, we develop a general approach to turn Monte Carlo samples with expectation $m$ into an unbiased estimate of $f(m)$. Specifically, we develop estimators that are based on randomly truncating the Taylor series expansion of $f$ and estimating the coefficients of the truncated series. We derive their properties and propose a strategy to set their tuning parameters – which depend on $m$ – automatically, with a view to make the whole approach simple to use. We develop our methods for the specific functions $f(x) = \log x$ and $f(x) = 1/x$, as they arise in several statistical applications such as maximum likelihood estimation of latent variable models and Bayesian inference for un-normalised models. Detailed numerical studies are performed for a range of applications to determine how competitive and reliable the proposed approach is.

## 1. INTRODUCTION

**1.1. General setting.** We consider the problem of unbiasedly estimating $f(m)$, for a smooth function $f$, using independent and identically distributed (IID) random variables $X_1, X_2, \dots$ with mean $m \in \mathbb{R}$. In particular, using the Taylor expansion of $f$ around a certain $x_0 \neq 0$,

$$ (1) \qquad f(m) = \sum_{k=0}^{\infty} \frac{f^{(k)}(x_0)}{k!} (m - x_0)^k = \sum_{k=0}^{\infty} \gamma_k \left( \frac{m}{x_0} - 1 \right)^k, \quad \gamma_k := \frac{f^{(k)}(x_0)}{k!} x_0^k, $$

we will employ the so-called ‘sum’ estimator (McLeish, 2011; Glynn and Rhee, 2014; Rhee and Glynn, 2015)

$$ (2) \qquad \hat{f} := \sum_{k=0}^{R} \frac{\gamma_k U_{R,k}}{\mathbb{P}(R \ge k)} = \sum_{k=0}^{\infty} \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)} \gamma_k U_{R,k}, $$

where $R$ is an integer-valued random variable,

$$ U_{R,k} = \sum_{r=k}^{\infty} U_{r,k} 1\{R=r\}, $$

and $U_{r,k}$ is an unbiased estimate of $(m/x_0 - 1)^k$, as clarified in the following assumption.

**Assumption 1.** Let $R$ be a random variable taking values in $\{0, 1, \dots\}$ such that $\mathbb{P}(R \ge k) > 0$ for all $k$. Let $\{U_{r,k} : k, r \in \mathbb{N}; k \le r\}$ be a collection of random variables, which are independent of $R$, and

$$ (3) \qquad \mathbb{E}[U_{r,k}] = \left(\frac{m}{x_0} - 1\right)^k. $$
Let $U_{0,0} = U_{1,0} = \cdots = 1$. Furthermore, we assume that $a_k := \sup_{r \ge k} \mathbb{E}[|U_{r,k}|] < \infty$ for $k > 0$, and $\sum_{k=0}^{\infty} |\gamma_k| a_k < \infty$.

These conditions ensure that $\hat{f}$ is an unbiased estimator, see the following lemma and Appendix A for a proof.

**Lemma 2.** Under Assumption 1, $\hat{f}$ is integrable and $\mathbb{E}[\hat{f}] = f(m)$ in (1).

A candidate for $U_{r,k}$ is the following construction of $U_{r,k}$ using $X_1, \dots, X_k$ only,

$$ (4) \qquad U_{r,k} = \prod_{i=1}^{k} (X_i/x_0 - 1), \quad 0 < k \le r. $$

Thus, $U_{r,k}$ unbiasedly estimates $(m/x_0 - 1)^k$ using the first $k$ univariates only. In Section 2.3, we propose a lower-variance estimator based on the whole sample $X_1, \dots, X_r$. Including the factor $x_0^k$ in the definition of $\gamma_k$ above will give simpler expressions for functions of interests, such as $f(m) = \log m$ and $f(m) = 1/m$. We shall also optimise this expansion point and henceforth, assume $x_0 \ne 0$.

## 1.2. Motivation and contributions.
We are particularly interested in the functions $f(x) = \log x$ and $f(x) = 1/x$. The former arises in applications such as maximum likelihood estimation where unbiased estimates of the likelihood are available, for example using importance sampling. Using the sum estimator framework, these can then be converted into unbiased estimates of the log-likelihood and its gradient. The function $f(x) = 1/x$ arises in Bayesian and frequentist estimation of models with un-normalised likelihoods, that is, when the likelihood involves an intractable normalising constant.

Various unbiased estimators have already been proposed in the literature for such functions; see Section 1.3 for a review. However, in our experience, their performance tends to be very sensitive to ‘tuning’ parameters such as the expansion point $x_0$ and the distribution of $R$. This is problematic for users, in particular in situations where one needs to routinely compute unbiased estimates of $f(m)$ at many different $m$ values. This is the case, for instance, in stochastic gradient methods, where at each iteration one needs to estimate the gradient of a different term; see Section 2.6 for details.

Thus, our motivation is to propose a general ‘turnkey’ approach to unbiased estimation which is both robust and easy to automate. Specifically, our contributions are as follows. We derive (Section 2) some results on the general properties of the randomly truncated Taylor-series based estimator, and use them to inform the choice of the truncation random variable $R$. We also discuss in detail the properties of the standard unbiased estimate (4) of $(m/x_0 - 1)^k$ and propose an alternative with potentially much smaller variance. In particular, our new estimate uses the whole collection of samples $X_1, \dots, X_r$ to construct the unbiased estimate $U_{r,k}$ of $(m/x_0 - 1)^k$ for $0 < k \le r$. We show the overall variance of the sum estimator will then decrease at a rate $O((\log \mathbb{E}[R])/\mathbb{E}[R])$. In contrast, the variance of the sum estimate with the standard choice (4) for the coefficient estimate does not diminish to zero as $\mathbb{E}[R]$ increases. This difference can be significant for some applications.

The second part of our study (Section 3) focuses on the specific functions $f(x) = \log x$ and $f(x) = 1/x$. For these functions, we discuss how to choose tuning parameters, $x_0$ and the distribution of $R$, to ensure that the resulting estimator has finite – and hopefully small – variance.
Section 4 presents numerical experiments that showcase the proposed approach in applications such as maximum likelihood estimation in latent variable models and Bayesian inference in un-normalised models. Section 5 concludes with a discussion of the merits and scope of the proposed approach. All the proofs appear in the Appendix.

1.3. **Related work.** Sum estimators and related debiasing techniques have gained traction since their introduction in McLeish (2011); Glynn and Rhee (2014); Rhee and Glynn (2015); for example, see Jacob and Thiéry (2015), Lyne et al. (2015), Agapiou et al. (2018) and Vihola (2018).

Some of these papers also consider single term estimators, which rely on a single randomly selected term to unbiasedly estimate the infinite series. In our particular context, single term estimators are less appealing since a single term of a Taylor expansion of $f(m)$ cannot give a good approximation. Therefore, we do not consider them further. We refer to, for example, Vihola (2018) for a more thorough discussion on the connection between single term estimators and sum estimators, and furthermore, their connection to multi-level Monte Carlo (Giles, 2008) in the context of stochastic differential equations.

The idea of applying a debiasing technique to a Taylor expansion appeared before in Blanchet et al. (2015), but without much discussion on how to make the corresponding estimators applicable in practical scenarios. Our work helps fill this gap.

Debiasing techniques specific to the log function have been developed in Luo et al. (2020) and Shi and Cornish (2021). However, as is shown in our numerical experiments, they generate estimates that may either have an infinite variance, or they have a CPU cost with infinite variance. In contrast, the CPU cost of our methods have a finite variance, and if properly tuned, the variance of the estimates are also finite.

The function $f(x) = \exp(x)$, which also has many important applications, is outside of the scope of this paper for reasons we explain later. We refer the reader to Papaspiliopoulos (2011) for a thorough treatment of debiasing techniques for this function.

1.4. **Notations and standing assumptions.** The random variables $X_1, X_2, \dots$ are independent and identically distributed (IID) with expected value $m$ and variance $\sigma^2$.

The expansion point $x_0 \in \mathbb{R} \setminus \{0\}$ in (1) is fixed initially, and thus so are the $\gamma_k$'s, but we discuss how to choose $x_0$, for a particular $m$, later on in Section 3.

Our results will be articulated when $R$ has the geometric distribution. We denote by Geometric($p$) the geometric distribution such that $\mathbb{P}(R = k) = (1-p)^k p$, $k = 0, 1, \dots$. Note that $\mathbb{E}[R] = (1-p)/p$ and $\mathbb{P}(R \ge k) = (1-p)^k$.

We say that an estimator has $O(R^k)$ complexity when its computation involves $O(r^k)$ operations, conditional on $R = r$. The estimators considered in this paper have all complexity at least $O(R)$, since they require generating $R$ random variables.

The notation $\hat{f} \in \mathcal{L}_2$ means $\mathbb{E}[|\hat{f}^2|] < \infty$, which implies that $\hat{f}$ has finite variance and is integrable, $\mathbb{E}[|\hat{f}||] < \infty$.
## 2. GENERAL PROPERTIES OF TAYLOR-BASED SUM ESTIMATORS

2.1. **Variance of the general sum estimator.** In this section, we consider the general estimator (2), without considering a particular choice for the $U_{r,k}$'s, beyond Assumption 1. Using this assumption,

$$ (5) \quad \mathbb{E}[\hat{f}|R=r] = \sum_{k=0}^{\infty} \gamma_k \frac{1\{r \ge k\}}{\mathbb{P}(R \ge k)} \mathbb{E}[U_{r,k}|R=r] = \sum_{k=0}^{\infty} \gamma_k \frac{1\{r \ge k\}}{\mathbb{P}(R \ge k)} \left(\frac{m}{x_0} - 1\right)^k. $$

See Appendix A.2 for the verification.

One may decompose its variance, assuming it is finite, through the law of total variance,

$$ (6) \quad \mathrm{var}[\hat{f}] = \mathrm{var}[\mathbb{E}[\hat{f}|R]] + \mathbb{E}[\mathrm{var}[\hat{f}|R]] $$

where the first term can be expressed as

$$ (7) \quad \mathrm{var}[\mathbb{E}[\hat{f}|R]] = \sum_{k=0}^{\infty} \gamma_k^2 \left(\frac{m}{x_0} - 1\right)^{2k} \left(\frac{1}{\mathbb{P}(R \ge k)} - 1\right) \\ \qquad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \gamma_k \gamma_l \left(\frac{m}{x_0} - 1\right)^{k+l} \left(\frac{1}{\mathbb{P}(R \ge k)} - 1\right). $$

See Appendix A.3 for the verification.

The term $\mathrm{var}[\mathbb{E}[\hat{f}|R]]$ has a simple interpretation, namely, it is the variance of $\hat{f}$ in the ideal case where the $X_i$'s and thus the $U_{r,k}$'s have zero variance, that is, $U_{r,k} = (m/x_0 - 1)^k$. Hence, this variance component is only due to the random truncation of the Taylor expansion.

To further analyse $\mathrm{var}[\mathbb{E}[\hat{f}|R]]$, we must make some assumptions on the $\gamma_k$'s.

**Assumption 3.** For the function $f: \mathbb{R} \to \mathbb{R}$, argument $m$ and expansion point $x_0$, the series expansion in (1) holds. Furthermore, (i) for some $c > 0$, $|\gamma_k| \le c$ for all $k = 0, 1, \dots$; (ii) $|\gamma_k| \exp(\alpha k) \to \infty$ for all $\alpha > 0$; (iii) $\beta_0 := |m/x_0 - 1| < 1$.

Note that the series (1) cannot converge if we impose (i) and (ii) without (iii).

**Proposition 4.** Under Assumptions 1 and 3, if $R \sim \text{Geometric}(p)$ with $p \in (0, 1 - \beta_0^2)$, then

$$ (8) \quad (f(m) - \gamma_0)^2 \left(\frac{p}{1-p}\right) \le \mathrm{var}[\mathbb{E}[\hat{f}|R]] \le \frac{c^2 \beta_0^2}{(1-\beta_0)^2} \times \frac{p}{1-p-\beta_0^2}. $$

See Appendix B.1 for a proof. Assumption 1 is needed to ensure integrability and thus that $\mathbb{E}[\hat{f}|R]$ is well defined. The condition $\beta_0 < 1$ of Assumption 3 and corresponding restriction on $p$ ensures the variance is finite.

The result above tells us that, under Assumption 3, one can take $R \sim \text{Geometric}(p)$ and $\mathbb{E}[\hat{f}|R]$ will have finite variance, provided $p < 1 - \beta_0^2$. In particular, the conditional variance behaves like $\mathcal{O}((\mathbb{E}[R])^{-1})$ since $\mathbb{E}[R] = 1/p - 1 \to \infty$ as $p \to 0$. We thus recover the standard Monte Carlo error rate, which is a variance that decreases as the inverse of the sample size.

**Remark 5.** We focus on the case $R \sim \text{Geometric}(p)$, since it is an easy-to-sample-from distribution; however, the upper bound established above holds for any $R$ such that $\mathbb{P}(R \ge k) \ge (1-p)^k$. Conversely, one can see by inspecting the proof that
a distribution for R with sub-geometric tails should typically lead to an infinite variance.

Of course, what remains to be seen is to what extent we can recover this ideal behaviour when we use the practical estimates, namely $U_{r,k}$'s, for $(m/x_0 - 1)^k$. This will be the topic of the next two sections.

**Remark 6.** Both our functions of interest fulfil Assumption 3: $\gamma_k = (-1)^{k-1}/k$ for $f(x) = \log x$, and $\gamma_k = (-1)^{k-1}$ for $f(x) = 1/x$. A function that does not fulfil Assumption 3 is $f(x) = \exp(x)$ as $\gamma_k$ decreases at a super-geometric rate. This makes it possible to consider for R some distribution with tails that are lighter than geometric, for example, the Poisson distribution, which is why we do not considered it further, but again see Papaspiliopoulos (2011) for a thorough treatment of this function.

Before moving on, we state the expression for the second term of decomposition (6), which we will study in the following sections. See Appendix A.4 for the derivation.

$$ (9) \quad \mathbb{E}[\mathrm{var}[\hat{f}|R]] = \sum_{k=0}^{\infty} \frac{\gamma_k^2}{\mathbb{P}(R \ge k)^2} \left\{ \sum_{r=k}^{\infty} \mathbb{P}(R=r) \mathrm{var}(U_{r,k}) \right\} \\ + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \left\{ \sum_{r=l}^{\infty} \mathbb{P}(R=r) \mathrm{cov}(U_{r,k}, U_{r,l}) \right\}. $$

2.2. Simple estimator. A standard choice (Papaspiliopoulos, 2011; Blanchet and Glynn, 2015) for $U_{r,k}$ is $U_{r,k}^S = U_{r,k}^S$, where

$$ U_{r,k}^S := \prod_{i=1}^{k} \left( \frac{X_i}{x_0} - 1 \right) $$

which actually does not depend on $r$. We denote by $\hat{f}^S$ the corresponding sum estimator.

Computing $\hat{f}^S$ requires generating $R$ random variables $X_i$, but only the first $k$ are used to compute the $k$-th term. It is easy to see that for any $0 \le k \le l$,

$$ \mathrm{cov}(U_{r,k}^S, U_{r,l}^S) = \left( \frac{\sigma^2}{x_0^2} + \left( \frac{m}{x_0} - 1 \right)^2 \right)^k \left( \frac{m}{x_0} - 1 \right)^{l-k} - \left( \frac{m}{x_0} - 1 \right)^{k+l}. $$

This leads to the following condition to ensure the finiteness of the variance of $\hat{f}^S$.

**Proposition 7.** In addition to Assumption 3, suppose $\beta < 1$ where

$$ (10) \qquad \beta^2 := \frac{\sigma^2}{x_0^2} + \left( \frac{m}{x_0} - 1 \right)^2 = \frac{\sigma^2}{x_0^2} + \beta_0^2, $$

and $R \sim \text{Geometric}(p)$ with $p \in (0, 1 - \beta^2)$. Then, $\hat{f}^S \in \mathcal{L}_2$ and

$$ (11) \qquad \mathbb{E}[\mathrm{var}(\hat{f}^S|R)] \le \frac{c^2(1+\beta)}{1-\beta} \times \frac{(1-p)}{1-p-\beta^2}. $$
See Appendix B.2 for a proof.

Note that, given $m$ and $\sigma^2$, it is always possible to find $x_0$ such that $\beta < 1$; in fact the minimiser of $\beta$ with respect to $x_0$ verifies this condition.

What is unsatisfactory about this result is that the right hand side does not vanish as $p \to 0$. This is not because the bound is not tight, as shown by the following example whose derivation is given in Appendix B.3. Our numerical experiments suggest that this is true for other functions of interest.

**Example 8.** Consider $f(x) = 1/x$, and thus $\gamma_k = (-1)^{k-1}$, and assume that we want to estimate $f(m)$ for $m > 0$. Then, we have

$$ \lim_{p \to 0} \mathbb{E} [\mathrm{var}[\hat{f}^S|R]] = \left(\frac{2}{m} - \frac{1}{x_0}\right) \left\{ \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} \right\}. $$

If $\sigma^2 > 0$, then $\beta^2 > \beta_0^2$, and since we need to select $x_0 > m/2$ to guarantee that the Taylor sum (1) converges, we have $\lim_{p \to 0} \mathbb{E}[\mathrm{var}[\hat{f}^S|R]] > 0$, showing that the variance of the simple estimator does not converge to 0 in general.

2.3. **Cycling estimator.** One obvious defect of the simple estimator is that computing $\hat{f}^S$ requires generating $R$ random variables, but $U_{R,k}^S$ estimates $(m/x_0 - 1)^k$ using only the first $k$ of these univariates. Utilising the full sample set $\{X_1, \dots, X_R\}$ to estimate $(m/x_0 - 1)^k$ could give more precise estimates.

We propose the following estimate obtained by averaging over circular shifts of the data, henceforth referred to as the ‘cycling’ estimate, (12)

$$ U_{r,k}^C := \frac{1}{r} \left[ \prod_{i=1}^{k} \left( \frac{X_i}{x_0} - 1 \right) + \prod_{i=2}^{k+1} \left( \frac{X_i}{x_0} - 1 \right) + \dots + \left( \frac{X_r}{x_0} - 1 \right) \prod_{i=1}^{k-1} \left( \frac{X_i}{x_0} - 1 \right) \right]. $$

Although it averages $r$ non-independent estimates of $(m/x_0 - 1)^k$, it should have lower variance than the simple estimator $U_{r,k}^S$. We denote by $\hat{f}^C$ the corresponding sum estimator.

The following proposition shows that the cycling estimate, indeed, typically has lower variance than the simple estimate.

**Proposition 9.** Assume the same conditions as Proposition 7, namely, Assumption 3, condition (10), and $R \sim \text{Geometric}(p)$ with $p \in (0, 1 - \beta^2)$. Then $\hat{f}^C \in \mathcal{L}_2$ and

$$ \mathbb{E} [\mathrm{var}(\hat{f}^C|R)] \le \frac{4c^2p\log(1/p)}{(1-p-\beta^2)^2} \left[ \beta^2 + 2 \frac{\beta_0(1-p)}{(1-\beta_0)^2} \right]. $$

See Appendix B.4 for a proof.

In other words, by taking $\mathbb{E}[R] \to +\infty$, the cycling estimate recovers, up to a log factor, the standard Monte Carlo rate, i.e. the variance is $\mathcal{O}((\log\{\mathbb{E}[R]\})/\mathbb{E}[R])$. Thus, both terms of the variance in (6) now decrease as the mean of $R$ increases.

A direct corollary, through Chebyshev's inequality, of the above result is that the cycling estimator converges in probability.

**Corollary 10.** Under the same assumptions as Proposition 9, $\hat{f}^C \xrightarrow{P} f(m)$ as $p \to 0$.
2.4. **Minimum variance unbiased estimator.** Instead of the cycling estimator, one could consider the following U-statistic:

$$ (13) \qquad U_{r,k}^{\text{MVUE}} := \frac{1}{\binom{r}{k}} \sum_{1 \le i_1 < \dots < i_k \le r} \left( \frac{X_{i_1}}{x_0} - 1 \right) \dots \left( \frac{X_{i_k}}{x_0} - 1 \right). $$

which has minimum variance among unbiased estimators of its expectation, see Hoeffding (1948). In spite of the combinatorial number of terms, it is still possible to compute the estimator in $O(r^2)$ operations, using a recursion, as pointed out by Sam Power in a personal communication, see Appendix D. However, repeating our variance analysis for the MVUE estimator appears challenging, indeed due to the combinatorial number of terms in (13), and is left for further study. For the moment, we present a numerical comparison in Appendix D, which happens to show cycling and MVUE performing similarly.

2.5. **Relative merits of the simple and cycling estimators.** The cycling estimator makes it possible to estimate $f(m)$ as accurately as required, and enjoys, up to a log factor, the same error rate as Monte Carlo estimates. On the other hand, it is more expensive to compute than the simple estimator. Although both estimators require generating $R$ random variables, computing the $U_{R,k}$'s ($0 \le k \le R$) for the cycling estimator requires $O(R^2)$ operations.

That said, in many practical scenarios, the $X_i$'s can be expensive to generate. For instance, in our third numerical experiment we run a sequential Monte Carlo (SMC) sampler to get a single $X_i$, which takes about 15 seconds. The CPU cost of the $O(R^2)$ operations for the cycling estimate is negligible compared to the CPU cost of generating the $X_i$'s, unless $R$ is extremely large. In this scenario, the cycling estimate is clearly preferable owing to the potentially large variance reductions achieved.

Another way to see the relative merits of both estimators is to look at their robustness relative to tuning parameters. This point will be discussed at length in Section 3. However, if our figure of merit is ‘variance times CPU cost’, which is a common figure of merit in Monte Carlo, then the cycling estimator $\hat{f}^C$ is clearly preferable since this figure of merit diverges at a log rate as $p \to 0$, whereas it diverges at a linear rate for the simple estimator.

2.6. **Unbiased estimators of gradients.** In this section, we further develop the simple and cycling estimators to unbiasedly estimate $\nabla_{\theta} f(m(\theta))$, the gradient of $f(m(\theta))$ when $m = m(\theta)$ depends on a parameter $\theta$. A typical application that utilises unbiased gradient estimates is maximum likelihood estimation; see Section 4 for an example.

Assuming $m(\theta)$ is differentiable, using (1),

$$ (14) \qquad \nabla_{\theta} f(m(\theta)) = \sum_{k=1}^{\infty} k \gamma_k \left( \frac{m(\theta)}{x_0} - 1 \right)^{k-1} \frac{\nabla_{\theta} m(\theta)}{x_0}, $$

and the corresponding ‘sum’ estimator, that is similar to (2) in construction, is

$$ (15) \qquad \widehat{\nabla_{\theta} f} := \sum_{k=1}^{\infty} \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)} W_{R,k} \frac{k \gamma_k}{x_0}, $$
where $R$ is a random variable taking values in $\{1, \dots\}$; $W_{R,k} = \sum_{r=k}^{\infty} W_{r,k} \mathbf{1}\{R=r\}$; $W_{r,k}$ is an unbiased estimate of $(m(\theta)/x_0 - 1)^{k-1}\nabla_{\theta}m(\theta)$ for $r \ge k$; and each and every $W_{r,k}$ is independent of $R$.

As was stated for $\hat{f}$ in Lemma 2, the gradient estimate can be shown to be integrable and unbiased when the estimates $W_{r,k}$ satisfy similar conditions of Assumption 1, that is, with univariates $U_{r,k}$ therein replaced with $W_{r,k}$, and $\gamma_k$ replaced with $k\gamma_k$. See Appendix A.1 for the proof of unbiasedness of $\hat{f}$ which can be easily adapted to the present gradient case.

We can define a simple estimator for (14), denoted by $\widehat{\nabla_{\theta}f}^S$, by replacing $W_{r,k}$ in (15) with

$$ (16) \qquad W_{r,k}^S := G_k \prod_{i=1}^{k-1} \left( \frac{X_i}{x_0} - 1 \right) $$

where $(X_1, G_1), (X_2, G_2), \dots, (X_r, G_r)$ is an IID collection of bivariates such that $\mathbb{E}[X_i] = m(\theta)$ and $\mathbb{E}[G_i] = \nabla_{\theta}m(\theta)$.

Instead of $W_{r,k}^S$, the cycling estimator, denoted by $\widehat{\nabla_{\theta}f}^C$, uses $W_{r,k}^C$ which is an average over ‘circular shifts’ of $W_{r,k}^S$,

$$ W_{r,k}^C := \frac{1}{r} \left[ G_k \prod_{i=1}^{k-1} \left( \frac{X_i}{x_0} - 1 \right) + G_{k+1} \prod_{i=2}^{k} \left( \frac{X_i}{x_0} - 1 \right) + \dots + G_{k-1} \prod_{i=1}^{k-2} \left( \frac{X_i}{x_0} - 1 \right) \right]. $$

As was done in the previous section, we can analyse the variance of both estimators using the law of total variance (6). To ensure that the second term of the variance decomposition, namely $\mathbb{E}[\mathrm{var}(\widehat{\nabla_{\theta}f}|R)]$, is finite we further need to assume that the second moment of $(X_1, G_1)$ is finite.

**Proposition 11.** At the point $\theta$, assume $(X_1, G_1)$ has finite second moments and let $\sigma^2(\theta) = \mathrm{var}(X_1)$. Furthermore, let Assumption 3 hold at $m = m(\theta)$; assume

$$ \beta^2 := \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 < 1; $$

and let $R \sim \text{Geometric}(p)$ with $p \in (0, 1 - \beta^2)$. Then both $\widehat{\nabla_{\theta}f}^S$ and $\widehat{\nabla_{\theta}f}^C$ are in $\mathcal{L}_2$. In addition,

$$ (1) \quad \mathrm{var}[\mathbb{E}[\widehat{\nabla_{\theta}f}|R]] = \mathcal{O}(p) \quad \text{as } p \to 0, \text{ for } \widehat{\nabla_{\theta}f} \in \{\widehat{\nabla_{\theta}f}^S, \widehat{\nabla_{\theta}f}^C\}. $$

$$ (2) \quad \text{As } p \to 0, \mathbb{E}[\mathrm{var}(\widehat{\nabla_{\theta}f}^S|R)] \text{ converges to a finite, possibly non-zero, value.} $$

$$ (3) \quad \mathbb{E}[\mathrm{var}(\widehat{\nabla_{\theta}f}^C|R)] = \mathcal{O}(p \log(1/p)) \quad \text{as } p \to 0. $$

The proof is in Appendix B.5.

The first statement implies the variance of the conditional expectation decreases as $\mathcal{O}(1/\mathbb{E}[R])$ as $\mathbb{E}[R] \to \infty$. The second statement confirms the variance of the simple estimate does not diminish to zero as $\mathbb{E}[R]$ increases; this is verified for $f(x) = \log(x)$. In contrast, the final statement confirms the cycle method’s variance diminishes at rate $\log(\mathbb{E}[R])/\mathbb{E}[R]$.
### 3. TUNING PARAMETERS

3.1. Tuning $x_0$. We now discuss how to set the tuning parameters $x_0$ and $p$ in practice. For the sake of simplicity, we focus on the cycling estimator, $\hat{f}^C$, and we assume $m > 0$, but adapting our discussion to either $\hat{f}^S$ or $m < 0$ should be straightforward.

Propositions 7 and 9 of Section 2 reveal the key considerations for selecting $x_0$ and parameter $p$ of $R \sim \text{Geometric}(p)$, appropriately. Firstly, to make sure that the variance of $\hat{f}^C$ is finite, $x_0$ must be chosen to ensure that $\beta^2 < 1$; see (10). Secondly, once $x_0$ is chosen, then $p$ has to be found from the interval $p \in (0, 1 - \beta^2)$. This order cannot be inverted in general because, if $p$ is too close to one, there may not be an $x_0$ that guarantees $\beta^2 < 1 - p$ as $\beta^2$ has a unique minimiser with respect to $x_0$.

We now go through how to make these critical choices. From our experience, this discussion is largely neglected in the literature. We will further illustrate in Section 4, with the aid of numerical examples, the importance of these tuning considerations, and the robustness of our tuning procedure below.

Recall, we focus the discussion for $m > 0$. We see that $\beta^2 < 1 \iff x_0 > x_0^{\min} := m/2 + \sigma^2/(2m)$. Rather than targeting $x_0^{\min}$, we target the $x_0$ that minimises $\beta^2$, namely

$$ (17) \qquad x_0^* := \arg \min_{x_0} \left\{ \frac{\sigma^2}{x_0^2} + \left( \frac{m}{x_0} - 1 \right)^2 \right\} = \frac{m^2 + \sigma^2}{m}. $$

Being the minimiser, $x_0^*$ automatically fulfils the condition for $\beta^2 < 1$. Incidentally, $x_0^* = 2x_0^{\min}$. But there are also other reasons why $x_0^*$ might be a good choice. The bound in Proposition 9 grows with $\beta^2$ and even diverges as $\beta^2 \to 1 - p$. So, it makes sense to minimise this quantity with respect to $x_0$. Also, minimising $\beta^2$ permits a larger $p$, and this could save on the computational cost to generate the estimate of $f(m)$.

The catch is that we do not know the values of $m$ and $\sigma^2$. The only general solution to this issue is to use a pilot run. That is, use an initial batch of $n_0$ random variables $X_i$ to obtain estimates $\hat{m}$ and $\hat{\sigma}^2$, plug these estimates into (17) to obtain $\hat{x}_0^*$ and use this estimate as the value for $x_0$.

Note however the following fundamental drawback of this approach. No matter how large the pilot run is, the probability that $\hat{x}_0^* < x_0^{\min}$, and equally $\beta^2 > 1$, will be non-zero, and the variance of $\hat{f}^C$ cannot be guaranteed to be finite.

To mitigate this issue, we need to find a way to make the probability $x_0 < x_0^{\min}$ as small as possible for the chosen $x_0$. We propose the following bootstrap (Efron and Tibshirani, 1993) approach. Generate a pilot sample of size $n_0$ and use it as follows:

i. Estimate $x_0^* = (m^2 + \sigma^2)/m$ as $\hat{x}_0^* = (\hat{m}^2 + \hat{\sigma}^2)/\hat{m}$ where $\hat{m}$ and $\hat{\sigma}^2$ are the empirical mean and variance of the pilot sample.

ii. Find the upper limit of a lower one-sided bootstrap $(1 - \alpha)$-confidence interval for $x_0^{\min}$, with $\alpha$ very small; note that $x_0^{\min} = x_0^*/2$.

iii. Set $x_0$ to the larger of these two values.

In this way, not only will the $P(x_0 > x_0^{\min})$ be close to one for the chosen $x_0$, the latter will potentially also be close to $x_0^*$ when the sample mean and variance are accurate.
In practice, we recommend using the percentile method to construct the bootstrap confidence interval, that is, the upper bound is the $(1 - \alpha)$-quantile of the bootstrap estimates of $x_0^{\min}$, as this method may work well even when the distribution of the $X_i$'s is not symmetric; for example, see Efron and Tibshirani (1993) for more background on bootstrap confidence intervals. In our experiments we set $\alpha = 0.01$. A second obvious recommendation is to favour a large $n_0$, as the computational budget permits.

**3.2. Setting p and the CPU budget.** Steps (i)-(iii) of Section 3.1 return an $x_0$ value such that $\mathbb{P}(x_0^{\min} < x_0) \ge 1 - \alpha$, where $x_0^{\min}$ is the plug-in estimate and $\mathbb{P}(\cdot)$ is computed with the empirical cumulative distribution function given by the pilot samples $\{X_i\}_{i=1}^{n_0}$. Let $\hat{\beta}^2$ be the corresponding plug-in estimate of $\beta^2$: in (10), use the found $x_0$ from Step (iii) and $(\hat{m}, \hat{\sigma}^2)$ from Step (i). Then – at this returned value for $x_0$ – evidently, we also have $\mathbb{P}(\hat{\beta}^2 < 1) \ge 1 - \alpha$. This motivates choosing $p$ as follows:

$$ \text{iv. Set } p = \min \left( 1 - \hat{\beta}^2, \frac{1}{n_0+1} \right). $$

The reason why we ‘align’ the choice of $p$ with the size of the pilot run $n_0$ is discussed next. Note that Geometric$(1/(n_0+1))$ has mean $n_0$. The resulting procedure is summarised in Algorithm 1 for the cycling estimator. The simple estimator is obtained by replacing $U_{r,k}^C$ with $U_{r,k}^S$.

1. Inputs: number of pilot runs $n_0$, a mechanism to sample $X_i$ with mean $m$, confidence level $\alpha$, $(\gamma_k)_{k \ge 0}$.

2. Tuning:

3. Sample $X_1, \dots, X_{n_0}$ and compute empirical mean $\hat{m}$ and variance $\hat{\sigma}^2$.

4. Estimate $\hat{x}_0^* = (\hat{m}^2 + \hat{\sigma}^2)/\hat{m}$.

5. Build a one-sided bootstrap $(1 - \alpha)$-confidence interval $(-\infty, u]$, for $\sigma^2/(2m) + m/2$.

6. Set $x_0 = \max\{u, \hat{x}_0^*\}$ and $\hat{\beta}^2 = \hat{\sigma}^2/x_0^2 + (\hat{m}/x_0 - 1)^2$.

7. Set $p = \min \left( 1 - \hat{\beta}^2, \frac{1}{n_0+1} \right)$.

8. Estimator:

9. Sample $R \sim \text{Geometric}(p)$ and, given $R = r$, sample $X_1, \dots, X_r$.

10. Compute $U_{r,k}^C$ as in (12).

11. Output: $\hat{f}^C = \sum_{k=0}^r \frac{\gamma_k U_{r,k}^C}{(1-p)^k}$.

**Algorithm 1:** Cycling Taylor-based sum estimator

Choosing $p$ amounts to choosing the number of variables $X_i$ that must be generated, and ultimately the CPU budget of $\hat{f}^C$ since, as we already mentioned, the generation of the $X_i$'s often dominates the CPU cost of $\hat{f}^C$. For both the cycling and simple estimators, we have seen in the previous section that a large pilot run of size $n_0$ is often required to choose $x_0$ well. It is thus illusory to try to use $\hat{f}^C$ in a low budget regime.

Put it in another way, the figure of merit we should now consider is the work-normalised variance WNV := $(n_0 + \mathbb{E}[R]) \times \text{var}[\hat{f}]$ (Glynn and Whitt, 1992). We do not consider the cost of building the estimators since, as discussed in Section 2.5,
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Comparison of the upper-bound of the work-normalised variance in a low, moderate and high variance scenario for $n_0 = 10$ as in the experiments in Section 4. The stars denote the value of $\mathbb{E}[R]$ minimising the upper-bound of the WNV and the vertical line $\mathbb{E}[R] = n_0 = 10$.</figcaption></figure>

we are mostly interested in scenarios in which the $X_i$'s are expensive to generate, see Section 4.3 for an example.

Figure 1 plots an upper bound on the work-normalised variance, obtained by replacing $\text{var}[\hat{f}]$ with the bounds given in Proposition 4 and 9 for the cycling estimator (Proposition 7 for the simple estimator) in the case $x_0 = x_0^*$, as a function of $\mathbb{E}[R]$. We consider $n_0 = 10$ since this is the value used in the experiments in Section 4. When $\mathbb{E}[R]$ gets too small, both bounds diverge at rate $O(1/(\mathbb{E}[R] - \beta^2))$; when $\mathbb{E}[R]$ tends to infinity, the bound for the cycling estimator grows slowly, at rate $O(\log \mathbb{E}[R])$, while that for the simple estimator diverges at rate $O(\mathbb{E}[R])$.

In practice, this means that the cycling estimator is more robust with respect to the choice of $p$: for larger values of $\mathbb{E}[R]$ the cost increases, but while the cycling estimator provides results with diminishing variance, the simple one fails to do so and the corresponding WNV explodes at a faster rate. This makes cycling preferable to the simple estimator.

We tried various ways to choose $p$ in an adaptive manner, by trying to minimise the upper bound on the work-normalised variance; see Appendix C for more details. However, we found in practice the following very basic strategy to work well. Set $\mathbb{E}[R] = n_0$, the size of the pilot run; see Algorithm 1. In this way, we fix the CPU budget beforehand to twice $n_0$ and given that $n_0 \gg 1$, this is often enough to ensure that we are in the right part of the curves in Figure 1. In the low and moderate variance cases, we also see that $\mathbb{E}[R] = n_0$ is not far from the value minimising the upper bound on WNV for the cycling and simple estimator.

## 4. NUMERICAL EXPERIMENTS

### 4.1. Toy latent variable model.
We consider the toy latent variable model (LVM) of Rainforth et al. (2018), also used in Shi and Cornish (2021). Conditioned on $\theta \in \mathbb{R}$, the observations $y_1, \dots, y_n$ are sampled independently as follows,

$$z_i|\theta \sim \mathcal{N}(\cdot; \theta 1_d, \text{Id}_d), \quad y_i|z_i \sim \mathcal{N}(\cdot; z_i, \text{Id}_d),$$

where $z_i$ and $y_i \in \mathbb{R}^d$.

We set $\theta = 1$ and sample a data set of size $n = 10^3$. The goal is to unbiasedly estimate the log-likelihood $\sum_{i=1}^n \log p(y_i|\theta)$ for $\theta = 1$, a quantity which we are able to compute exactly for this model.
To do so, we set $f(m) = \log m$ and consider each of the $n$ data points separately. For a given $y \in \{y_1, \dots, y_n\}$, we apply our approach as follows: conditional on $R = r$, we sample $Z_1, \dots, Z_r \sim p(z|\theta)$ and set $X_j = p(y|Z_j, \theta)$ for $j = 1, \dots, r$, which is evidently an unbiased estimate of $p(y|\theta) = \int p(y|z, \theta)p(z|\theta)dz$. In this particular case, we can exactly compute the moments $m = \mathbb{E}[X_j] = p(y|\theta)$ and $\sigma^2 = \text{var}[X_j]$; the optimal expansion point $x_0^* = m + \sigma^2/m$; and the corresponding value for $\beta^2$. See Appendix E.1 for their expressions.

We use this toy model to compare the performance of the simple and cycling estimators of $\log p(y|\theta)$ with that of the multilevel Monte Carlo (MLMC) estimator of Shi and Cornish (2021). They use MLMC to modify the SUMO approach of Luo et al. (2020) to return an unbiased estimate of $\log p(y|\theta)$ with finite variance. Thus, we do not consider the SUMO estimator here since its variance is infinite. Furthermore, Shi and Cornish (2021, Section 7.1) show, empirically, that MLMC outperforms SUMO.

The estimator proposed by Shi and Cornish (2021) employs MLMC (Blanchet et al., 2019) to debias the estimator of $\log p(y|\theta)$ given by Importance Weighted Autoencoders (IWAE; Burda et al. (2016)). The MLMC estimator is

$$ (18) \qquad \widehat{\log p(y|\theta)}^{\text{MLMC}} = I_j + \sum_{k=0}^{\infty} \frac{\Delta_k \mathbb{1}\{\tilde{R} \ge k\}}{\mathbb{P}(\tilde{R} \ge k)}, $$

where $\tilde{R} \sim \text{Geometric}(0.6)$, $I_j, j \ge 0$, is the IWAE estimate with proposal $q(\cdot; y)$ and $2^j$ samples:

$$ I_j := \log \frac{1}{2^j} \sum_{h=1}^{2^j} w(Z_h; y, \theta), \quad \text{with } w(Z_h; y, \theta) := \frac{p(Z_h, y|\theta)}{q(Z_h; y)}, \quad Z_1, \dots, Z_{2^j} \sim q(\cdot; y). $$

and $\Delta_k := I_{j+k+1} - (I_{j+k}^E + I_{j+k}^O)/2$; with $I_{j+k}^E$ (resp. $I_{j+k}^O$) being the IWAE estimate that only uses the even-indexed (resp. odd-indexed) samples. We consider the same proposal as in Shi and Cornish (2021), i.e. $\mathcal{N}(\cdot; (y + \theta^* 1_d)/2, \frac{2}{3} \mathbb{1}_{|d|})$, where $\theta^*$ denotes the maximum likelihood estimator for $\theta$, which is shown to be optimal within a certain family in Rainforth et al. (2018).

The expected CPU cost, i.e. the number of generated $X_j$, is $\mathbb{E}[R] = 1/p - 1$ for our Taylor-based estimators, and $\mathbb{E}[2^{j+1+\tilde{R}}]$ for the MLMC approach of Shi and Cornish (2021). We adjust $p$ and $j$ so that both expectations are equal to a certain pre-specified budget $C$. We consider two values for $C$, namely $C = 6$ and $C = 96$.

Figure 2 shows 100 estimates of each type (simple, cycling, MLMC) for $d = 2$ (left column) and $d = 5$ (right column) against their computational cost, for an expected cost $C$ of 6 (top line) and 96 (bottom line) samples per data point. The actual computational cost (reported on the y-axis) is the total number of variables $Z_j$ generated in order to estimate all the terms $\log p(y_i|\theta)$ for $i = 1, \dots, n$.

This experiment shows that the cost of MLMC has large empirical variance. In fact, the true variance of the cost is infinite, see Appendix E.1 for a proof of this statement. In contrast, the cost of the simple and cycling estimators remains close to its expectation.

The simple and cycling estimators outperform MLMC in terms of accuracy, with the simple estimator having larger variance than the cycling estimator for higher computational cost. This is particularly evident in the first plot of the second row, in which we compare the estimates for $\mathbb{E}[R] = 96$, in this case the cloud of
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Estimates of the marginal log-likelihood for a full data set of size $n = 10^3$ plotted against computational cost (measured in number of samples $Z_j$ drawn) for an expected cost $C$ of 6 and 96 samples per data point. The first column shows the results for $d = 2$ while the second one for $d = 5$. The vertical dashed line denotes the expected cost while the horizontal one denotes the true value of $\sum_{i=1}^{n} \log p(y_i|\theta)$.</figcaption></figure>

estimates obtained with cycling is considerably more concentrated around the true value (the horizontal line). This is consistent with the theoretical results established in Section 2 which show that the variance of the simple estimator does not decay to 0 when the computational cost increases.

As a note of caution, our Taylor-based estimators may give very poor results for a larger dimension $d$; see Appendix E.1 for extra results when $d = 20$. Simple calculations show that, in that case, one would need to take $E[R] \gg 10^8$ to ensure that $p < 1 - \beta^2$ and that the estimators have finite variance. In this case, the inputs have very large variance, and presumably the only way to fix this issue is to come up with a different strategy to unbiasedly estimate $p(y|\theta)$, so that the corresponding estimates have much lower variance.

### 4.2. Independent component analysis.
We consider a LVM associated with probabilistic independent component analysis (ICA). The data $y_1, \dots, y_n$ are $d_y$-dimensional vectors obtained from

$$ (19) \qquad y = \sum_{j=1}^{d_z} z_j a_j + \sigma \varepsilon, $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Estimated decomposition images (parameter vector $A$ of $\theta$) with SAEM, MLMC-SGD and Taylor-based sum estimators for the censored logistic ICA model. The first column shows the true value of $a_1, a_2$.</figcaption></figure>

where $z_j$ are the latent random variables; $a_j$'s are parameter vectors of dimension $d_y$; $\varepsilon$ is a $d_y$-dimensional standard Gaussian vector and $\sigma > 0$. Furthermore, $z_j = b_j\zeta_j$ where $b_j \sim \text{Bernoulli}(\alpha)$ and $\zeta_j \sim \text{Logistic}(1/2)$ for some fixed $\alpha$. The aim is to estimate the parameter vector $\theta = (A, \sigma)$, with $A := (a_1, \dots, a_{d_z})$, via maximum marginal-likelihood estimation, that is, we want to find the $\theta$ that maximises $\sum_{i=1}^n \log p(y_i|\theta)$.

We repeat the first experiment in Allassonnière and Younes (2012). We generate $n = 100$ data points from (19) with the following model choices: $d_z = 2$, $\sigma = 0.5$, $\alpha = 0.8$. Vectors $a_1$ and $a_2$ are the 16 × 16 images in the first panel of Figure 3. Note that $d_y = 256$ and the parameter space has dimension $d_\theta = 2 \times 256 + 1$.

To estimate $\theta$, Allassonnière and Younes (2012) use a stochastic approximation expectation maximisation (SAEM) algorithm. We will instead use stochastic gradient descent (Robbins and Monro, 1951), implemented with our unbiased estimators of $\nabla_\theta f(m(\theta)) = \nabla_\theta \log p(\theta|y)$. We call our two implementations simple-SGD and cycling-SGD. We also implement SGD with the gradient unbiasedly estimated using the MLMC strategy of Shi and Cornish (2021) (MLMC-SGD). The bivariates $(X_i, G_i)$ (for $i = 1, 2, \dots$) of Section 2.6 are generated using importance sampling; and the same importance sampling proposal is used for MLMC. For MLMC, we adjust its parameters $p$ and $j$, defined in (18) with $\tilde{R} \sim \text{Geometric}(p)$, so that $\mathbb{E}[2^{j+1+\tilde{R}}] = n_0$. Similarly, $p$ for the Geometric random truncation of the simple and cycling gradient estimates is found by setting $\mathbb{E}[R] = n_0$. See Appendix E.3 for all the details.

Table 1 shows the mean square error (MSE) of the estimates and their cost over 100 repetitions. Figure 3 shows the estimated images (parameter $A$ of $\theta$) from one run of each algorithm. SAEM is more costly than the other methods since it scans the whole data set of size $n = 100$ at each iteration, while the other SGD methods use randomly selected mini-batches of size 2. The overhead caused by cycling is small compared to the cost of running SGD for 5000 steps. Although the final MSE is similar for the four algorithms, on a visual inspection (Figure 3), SGD seems to yield better $A$ estimates than SAEM.
<table><caption>TABLE 1. Comparison of reconstruction accuracy and cost for the censored logistic ICA model. The results are averaged over 100 repetitions.</caption><thead><tr><th>Method</th><th>MSE(σ)</th><th>MSE(A)</th><th>E [Nsamples]</th><th>var (Nsamples)</th><th>runtime (s)</th></tr></thead><tbody><tr><td>SAEM</td><td>7.00 · 10<sup>-7</sup></td><td>0.22</td><td>—</td><td>—</td><td>33</td></tr><tr><td>MLMC</td><td>2.95 · 10<sup>-6</sup></td><td>0.29</td><td>12</td><td>19311</td><td>11</td></tr><tr><td>simple</td><td>2.97 · 10<sup>-6</sup></td><td>0.27</td><td>10</td><td>90</td><td>10</td></tr><tr><td>cycling</td><td>6.80 · 10<sup>-5</sup></td><td>0.24</td><td>10</td><td>90</td><td>14</td></tr></tbody></table>

We also observe that there seems to be a small difference between the results obtained with the simple and the cycling estimator. Both are unbiased, however, while running our experiments we noticed that the higher variance of the simple estimator requires smaller learning rates to counterbalance the gradient estimate.

### 4.3. Exponential random graphs.
Exponential random graph (ERG) models are one of the classical examples of models with a doubly intractable posterior distributions. For this model, the data is the observed edge connectivity of an n-node undirected graph. In particular, one observes $y = (y_{ij})_{i<j}$, where $y_{ij} = 1$ (resp. 0) if nodes i and j are connected (resp. are not). The likelihood is, given a parameter $\theta$ and canonical statistics $s(y)$, both in $\mathbb{R}^k$,

$$ (20) \qquad p(y|\theta) = \exp\{\theta^T s(y)\} / \mathcal{Z}(\theta). $$

The normalising constant $\mathcal{Z}(\theta)$ is a sum over $2^{\binom{n}{2}}$ terms, and is therefore intractable.

We construct an unbiased estimate of $1/\mathcal{Z}(\theta)$, for any $\theta$, which may then be combined with the remaining tractable terms of (20) to perform Bayesian inference and model choice. We consider the same example as in Caimo and Friel (2011): an $n = 16$ node graph; the statistic is $s(y) = (\sum_{i<j} y_{ij}, \sum_{i<j<k} y_{ij}y_{jk})$, which counts the number of edges and of two stars; $\theta \in \mathbb{R}^2$ is assigned the prior $N_2(0_2, 30I_2)$. The data – the observed connectivity – is the Florentine family business graph (Padgett and Ansell, 1993).

For a given $\theta$, we will unbiasedly estimate $1/\mathcal{Z}(\theta)$ using the simple and cycling estimator applied to the function $f(m) = 1/m$. The unbiased estimates of $m = \mathcal{Z}(\theta)$, or inputs $X_i$ used in the simple and cycling methods, are obtained by running the waste-free variant (Dau and Chopin, 2022) of tempered SMC; see Appendix E.4 for more details on this algorithm. In this example, the cost of generating the inputs with SMC is significantly larger than the cost of assembling the cycling estimate.

We first fix $\theta = (-3.15, 0.58)$, a value in the high probability region of the posterior, and compare the performance of the simple and the cycling estimator.

We consider two scenarios, one in which the unbiased estimates of $\mathcal{Z}(\theta)$ have moderate variance and one in which the variance is low. See Appendix E.4 for details on the moderate and low variance regimes. Basically, the latter is obtained by taking 10 times more particles in SMC, which leads to a (roughly) ten times smaller variance for the inputs.

We select the value of $x_0$ using the strategy in Section 3.1 with $n_0 = 10$. Once $x_0$ has been selected, we identify the corresponding $\hat{\beta}^2$ as discussed in Section 3.2, see Algorithm 1 for a summary. This gives us a range for $p \in (0, 1 - \hat{\beta}^2)$. In the moderate variance regime we obtain $p \in (0, 0.7)$ while for the low variance regime
<table><caption>TABLE 2. Work-normalised variance and proportion of negative estimates out of $10^3$ replicates of $\tilde{Z}(\theta)^{-1}$ for the simple and the cycling estimator when the variance of the estimates of $Z(\theta)$ is moderate and low.</caption><thead><tr><th rowspan="3"></th><th colspan="4">Simple</th><th colspan="4">Cycling</th></tr><tr><th colspan="2">Moderate</th><th colspan="2">Low</th><th colspan="2">Moderate</th><th colspan="2">Low</th></tr><tr><th>WNV</th><th>P(-)</th><th>WNV</th><th>P(-)</th><th>WNV</th><th>P(-)</th><th>WNV</th><th>P(-)</th></tr></thead><tbody><tr><th>E[R] = 1</th><td>111</td><td>0.026</td><td>106</td><td>0.006</td><td>102</td><td>0.028</td><td>77</td><td>0.002</td></tr><tr><th>E[R] = 10</th><td>61</td><td>0.027</td><td>83</td><td>0.001</td><td>24</td><td>0.01</td><td>26</td><td>0.0</td></tr></tbody></table>

we have $p \in (0, 0.95)$. In both cases, it is sufficient to set $E[R] = 1$ (i.e. $p = 0.5$) to ensure that both the simple and the cycling estimator have finite variance.

Since the cost of obtaining the estimates of $Z(\theta)$ is significantly larger than that of computing the estimates, selecting $E[R] = 1$ might be wasteful as our theoretical results in Section 2 guarantee that the cycling estimator has a lower variance than the simple estimator for large values of $E[R]$. Thus, we also consider $E[R] = 10 = n_0$ in our simulations as recommended in Section 3.2.

Table 2 compares the simple and cycling estimators in these four different scenarios (two variance regimes, two values of $E[R]$) according to two metrics estimated on $10^3$ independent runs: work-normalised variance (relative variance times CPU cost), and proportion of negative $1/Z(\theta)$ estimates. The results show that (a) the cycling estimator systematically outperforms the simple one; and (b) increasing $E[R]$ does lead to better performance for the cycling estimator but not for the simpler one, as expected. The CPU cost of cycling is less than 0.01 seconds higher than that of the simple estimator.

To recover the posterior distribution and the marginal likelihood, we use importance sampling, with a Gaussian proposal $q$, from which we generate $n = 1024$ vectors $\theta_i$ through randomised quasi-Monte Carlo. Using quasi-Monte Carlo ensures that the variability of our importance sampling estimates is dominated by the randomness of the estimates of $1/Z(\theta)$. See the appendix for more details on the proposal and on randomised quasi-Monte Carlo.

The importance sampling weights are

$$ (21) \qquad w_i = \frac{p(\theta_i) \exp\{\theta_i^T s(y)\}}{q(\theta_i)} \times \frac{1}{\widehat{Z}(\theta_i)}, \quad \theta_i \sim q. $$

The average of these noisy weights gives us an unbiased estimate of the marginal likelihood, $p(y) = \int p(\theta)p(y|\theta)d\theta$. To the best of our knowledge, this is the first study that uses unbiased estimates of the marginal likelihood of an ERG model.

For each $\theta_i$, we use exactly the same approach to estimate $1/Z(\theta_i)$ as the one which led to best performance in the previous section, namely, $n_0 = 10$, $E[R] = 10$, and waste-free SMC run in the “moderate variance” regime.

Figure 4 shows the histograms obtained with the importance sampling approximations of the posterior distribution $p(\theta|y)$; ‘simple’ and ‘cycling’ correspond to (21) with $1/Z(\theta)$ estimated with the simple and cycling methods respectively. The simple (resp. cycling) estimator returns 35 (resp. 10) negative weights out of 1024 samples. The effective sample size (ESS), i.e. $(\sum_i w_i)^2 / \sum_i (w_i)^2$ is 6% (resp. 31%) when using the simple (resp. cycling) estimator. Using the cycling estimator here,
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Bivariate weighted histograms approximating the posterior distributions obtained with the simple and the cycling estimator using $n = 1024$ samples from proposal $q$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Distribution of ESS and evidence (marginal likelihood) for 100 repetitions of the posterior approximation using the simple and the cycling estimator.</figcaption></figure>

approximately, yields a fivefold increase in computational efficiency with little additional cost.

To assess the variability of our approach we repeat the estimation of $1/Z(\theta)$, for each $\theta_i$, one hundred times, resulting in 100 different posterior approximations obtained from the same set of $n = 1024$ samples obtained from the proposal. For each repetition, we compute the ESS, and the unbiased approximation of the model evidence $p(y)$ given by $n^{-1} \sum_{i=1}^{n} w(\theta_i)$ (Figure 5). The ESS obtained with the cycling estimator is generally higher than that obtained with the simple estimator, and is above 0.2 in 74 out of the 100 repetitions. The average number of negative weights for a repetition is 44 for the simple estimator and 18 for the cycling estimator.

Low values of the ESS correspond to outliers in the distribution of the evidence, in particular, negative values of $p(y)$ are more frequent with the simple estimator than with the cycling estimator. The resulting variance for the estimates of $p(y)$ is 6 times larger for the simple estimator than for the cycling estimator.

## 5. CONCLUSION

We develop a general strategy for unbiasedly estimating $f(m)$, for a smooth function $f$, from random variables $X_i$ with expectation $m$. Inspired by the simple estimator (Section 2.2) – which is widely used in the literature – we proposed a new cycling estimator (Section 2.3), and obtained conditions that guarantee finite
estimator and execution time (CPU cost) variance for both methods. This is an improvement over previous approaches. In particular, the multi-level Monte Carlo method of Luo et al. (2020) returns estimates with infinite variance, while the one of Shi and Cornish (2021) has a running time with infinite variance. The latter issue may be particularly problematic in a parallel implementation since the completion time is the maximum of all threads being run.

As pointed out by a referee, the fact that the execution cost is random may be a drawback in certain implementation scenarios, even if it has finite variance. In such a case, one may prefer using a biased estimate whose execution cost is deterministic. Such an approach has less theoretical support, but seems to work well in certain practical applications, e.g., to approximate gradients in variational auto-encoders (Burda et al., 2016).

We show that the cycling estimator outperforms the simple one both theoretically, see Proposition 9 and Example 8 and empirically, see Section 4, with potentially only a little more computational cost. We also discussed the sensitivity of these methods to their tuning parameters, an issue we found to be largely overlooked in the literature. We gave precise guidelines for selecting the expansion point $x_0$ and the distribution of the random variable $R$ that determines the random truncation (Sections 3.1 and 3.2).

Our experiments show that, in order to build reliable unbiased estimates of $f(m)$ or its gradient, one needs input variables $X_i$ with low or moderate variance. In fact, the minimal $\beta^2$ in Proposition 7, achieved for $x_0^*$ in (17), is $(\beta^2)^* = \sigma^2/(\sigma^2 + m^2)$. Note that $1 - (\beta^2)^*$ is indeed the inverse of relative variance of the $X_i$'s. As $\sigma^2/m^2 \to \infty$, we need to select $p \in (0, 1 - (\beta^2)^*) \to \{0\}$, hence requiring an infinite computational cost. In this case, we recommend dedicating the computational effort towards building lower variance estimates of $m$.

Even if the $X_i$'s have low variance, a poorly chosen $x_0$ can cause the estimators of $f(m)$ to have high variance. While in low dimensional problems one can get away with a poorly tuned $x_0$, we found that for high dimensional applications, e.g., the ICA example in Section 4.2, it is crucial that the input variables $X_i$'s have low variance and to fine tune $x_0$, as recommended in Section 3.1. Once $x_0$ has been selected, we describe a simple strategy to select $p$ to guarantee that $E[R]$ is sufficiently large to have finite and small variance. This strategy works well for the experiments we considered, however, in some cases it might be possible to build optimal distributions for particular $f$'s (see, e.g. Rhee and Glynn (2015); Cui et al. (2021); Zheng et al. (2023)).

As a by-product of our study, we provide a method to unbiasedly estimate the model evidence in doubly intractable models. Previous attempts in the context of ERG models have focused on pseudo-likelihood approaches in which the intractable likelihood is replaced by a tractable approximation (Bouranis et al., 2018). Also, using a combination of path sampling and kernel density estimation, which provides biased estimators and are only suitable in low-dimensional settings (Caimo and Friel, 2013, Section 4).

## ACKNOWLEDGMENTS

Research of the first two authors is partly supported by the French Agence Nationale de la Recherche (ANR) through grant EPISVACAGE. SSS holds the Tibra Foundation professorial chair and gratefully acknowledges research funding
as follows: This material is based upon work supported by the Air Force Office of Scientific Research under award number FA2386-23-1-4100. The authors would like to thank Otmane Sakhi, Adrien Corenflos Sam Power and two referees for helpful comments at various stages of this project.

## REFERENCES

Agapiou, S., Roberts, G. O., and Vollmer, S. J. (2018). Unbiased Monte Carlo: Posterior estimation for intractable/infinite-dimensional models. *Bernoulli*, 24(3):1726–1786.

Allassonnière, S. and Younes, L. (2012). A stochastic algorithm for probabilistic independent component analysis. *Annals of Applied Statistics*, 6(1):125–160.

Beatson, A. and Adams, R. P. (2019). Efficient optimization of loops and limits with randomized telescoping sums. In *International Conference on Machine Learning*, pages 534–543. PMLR.

Billingsley, P. (1995). *Probability and Measure*. John Wiley & Sons.

Blanchet, J. H., Chen, N., and Glynn, P. W. (2015). Unbiased Monte Carlo computation of smooth functions of expectations via Taylor expansions. In *2015 Winter Simulation Conference (WSC)*, pages 360–367. IEEE.

Blanchet, J. H. and Glynn, P. W. (2015). Unbiased Monte Carlo for optimization and functions of expectations via multi-level randomization. In *2015 Winter Simulation Conference (WSC)*, pages 3656–3667. IEEE.

Blanchet, J. H., Glynn, P. W., and Pei, Y. (2019). Unbiased multilevel Monte Carlo: Stochastic optimization, steady-state simulation, quantiles, and other applications. *arXiv preprint arXiv:1904.09929*.

Boas Jr, R. P. (1965). Tannery's Theorem. *Mathematics Magazine*, 38(2):66–66.

Bouranis, L., Friel, N., and Maire, F. (2018). Bayesian model selection for exponential random graph models via adjusted pseudolikelihoods. *Journal of Computational and Graphical Statistics*, 27(3):516–528.

Burda, Y., Grosse, R., and Salakhutdinov, R. (2016). Importance weighted autoencoders. In *4th International Conference on Learning Representations (ICLR)*.

Caimo, A. and Friel, N. (2011). Bayesian inference for exponential random graph models. *Social Networks*, 33(1):41–55.

Caimo, A. and Friel, N. (2013). Bayesian model selection for exponential random graph models. *Social Networks*, 35(1):11–24.

Cui, Z., Lee, C., Zhu, L., and Zhu, Y. (2021). On the optimal design of the randomized unbiased Monte Carlo estimators. *Operations Research Letters*, 49(4):477–484.

Dau, H.-D. and Chopin, N. (2022). Waste-free sequential Monte Carlo. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 84(1):114–148.

Delyon, B., Lavielle, M., and Moulines, E. (1999). Convergence of a stochastic approximation version of the EM algorithm. *Annals of Statistics*, pages 94–128.

Efron, B. and Tibshirani, R. J. (1993). *An introduction to the bootstrap*, volume 57 of *Monographs on Statistics and Applied Probability*. New York, NY: Chapman & Hall.

Giles, M. B. (2008). Multilevel Monte Carlo path simulation. *Operations Research*, 56(3):607–617.

Glynn, P. W. and Rhee, C.-h. (2014). Exact estimation for Markov chain equilibrium expectations. *Journal of Applied Probability*, 51(A):377–389.
Glynn, P. W. and Whitt, W. (1992). The asymptotic efficiency of simulation estimators. *Operations Research*, 40(3):505-520.

Hoeffding, W. (1948). A non-parametric test of independence. *Annals of Mathematical Statistics*, 19(4):546-557.

Jacob, P. E. and Thiéry, A. H. (2015). On nonnegative unbiased estimators. *Annals of Statistics*.

Kuhn, E., Matias, C., and Rebafka, T. (2020). Properties of the stochastic approximation EM algorithm with mini-batch sampling. *Statistics and Computing*, 30(6):1725-1739.

Luo, Y., Beatson, A., Norouzi, M., Zhu, J., Duvenaud, D., Adams, R. P., and Chen, R. T. Q. (2020). SUMO: Unbiased estimation of log marginal probability for latent variable models. In *International Conference on Learning Representations (ICLR)*.

Lyne, A.-M., Girolami, M., Atchadé, Y., Strathmann, H., and Simpson, D. (2015). On Russian roulette estimates for Bayesian inference with doubly-intractable likelihoods. *Statistical Science*, 30(4):443-467.

McLeish, D. (2011). A general method for debiasing a Monte Carlo estimator. *Monte Carlo Methods and Applications*, 17(4):301-315.

Owen, A. B. (1997). Scrambled net variance for integrals of smooth functions. *Annals of Statistics*, 25(4):1541-1562.

Padgett, J. F. and Ansell, C. K. (1993). Robust action and the rise of the Medici, 1400-1434. *American Journal of Sociology*, 98(6):1259-1319.

Papaspiliopoulos, O. (2011). Monte Carlo probabilistic inference for diffusion processes: A methodological framework. In Barber, D., Cemgil, A. T., and Chiappa, S., editors, *Bayesian Time Series Models*, pages 82–103. Cambridge University Press.

Polyak, B. T. and Juditsky, A. B. (1992). Acceleration of stochastic approximation by averaging. *SIAM Journal on Control and Optimization*, 30(4):838-855.

Rainforth, T., Kosiorek, A., Le, T. A., Maddison, C., Igl, M., Wood, F., and Teh, Y. W. (2018). Tighter variational bounds are not necessarily better. In *International Conference on Machine Learning*, pages 4277–4285. PMLR.

Rhee, C.-h. and Glynn, P. W. (2015). Unbiased estimation with square root convergence for SDE models. *Operations Research*, 63(5):1026-1043.

Robbins, H. and Monro, S. (1951). A stochastic approximation method. *Annals of Mathematical Statistics*, pages 400-407.

Shi, Y. and Cornish, R. (2021). On multilevel Monte Carlo unbiased gradient estimation for deep latent variable models. In *International Conference on Artificial Intelligence and Statistics*, pages 3925-3933. PMLR.

Vihola, M. (2018). Unbiased estimators and multilevel Monte Carlo. *Operations Research*, 66(2):448-462.

Zheng, C., Pan, J., and Wang, Q. (2023). Optimal distributions for randomized unbiased estimators with an infinite horizon and an adaptive algorithm. preprint arXiv:2304.07797.
APPENDIX A. INTEGRABILITY AND MOMENT EXPRESSIONS

A.1. **Proof of Lemma 2.** We will show

$$ (22) \qquad \mathbb{E} \left[ \sum_{k=0}^{\infty} \frac{1\{R \geq k\}}{\mathbb{P}(R \geq k)} |\gamma_k| |U_{R,k}| \right] < \infty, $$

which implies the infinite series on the right-hand side of (2) is, almost surely, absolutely summable. Thus, $\hat{f}$ exists almost surely and is integrable. To show it has the desired mean, as (22) holds, we can apply Fubini's Theorem for general functions (Billingsley, 1995, Theorem 18.3) to swap the expectation and the sum and obtain

$$ \mathbb{E}[\hat{f}] = \sum_{k=0}^{\infty} \frac{\mathbb{E}[1\{R \geq k\}U_{R,k}]}{\mathbb{P}(R \geq k)} \gamma_k = \sum_{k=0}^{\infty} \gamma_k \left( \frac{m}{x_0} - 1 \right)^k. $$

To show (22), we note that all the terms in the sum are positive. Fubini's Theorem for non-negative functions (Billingsley, 1995, Theorem 18.3) allows the swap of the expectation and the sum,

$$ \mathbb{E} \left[ \sum_{k=0}^{\infty} \frac{1\{R \geq k\}}{\mathbb{P}(R \geq k)} |\gamma_k| |U_{R,k}| \right] = \sum_{k=0}^{\infty} \frac{\mathbb{E}[1\{R \geq k\}|U_{R,k}|]}{\mathbb{P}(R \geq k)} |\gamma_k| \leq \sum_{k=0}^{\infty} |\gamma_k| a_k < \infty, $$

where we used

$$ \begin{aligned} \mathbb{E}[1\{R \geq k\}|U_{R,k}|] &= \sum_{r=k}^{\infty} \mathbb{E}[1\{R=r\}|U_{r,k}|] \\ &= \sum_{r=k}^{\infty} \mathbb{P}(R=r)\mathbb{E}[|U_{r,k}|] \\ &\leq a_k \sum_{r=k}^{\infty} \mathbb{P}(R=r). \end{aligned} $$

Note that the independence of $R$ and $U_{r,k}$ was invoked in the second line.

A.2. **Verification of (5).** The conditional expectation exists, and is integrable, since $\hat{f}$ itself is integrable. By definition, $\mathbb{E}[\hat{f}|R=r] = \mathbb{E}[\hat{f}1\{R=r\}]/\mathbb{P}(R=r)$, provided $\mathbb{P}(R=r)>0$. It follows from (22) that

$$ (23) \qquad \mathbb{E} \left[ \sum_{k=0}^{\infty} \frac{1\{R \geq k\}}{\mathbb{P}(R \geq k)} |\gamma_k| |U_{R,k}| 1\{R=r\} \right] < \infty. $$

Since (23) holds, we can apply Fubini's Theorem for general functions (Billingsley, 1995, Theorem 18.3) and obtain

$$ \mathbb{E}[\hat{f}|R=r] = \sum_{k=0}^{\infty} \gamma_k \frac{1\{r \geq k\}}{\mathbb{P}(R \geq k)} \mathbb{E}[U_{r,k}|R=r] = \sum_{k=0}^{\infty} \gamma_k \frac{1\{r \geq k\}}{\mathbb{P}(R \geq k)} \left( \frac{m}{x_0} - 1 \right)^k. $$

Note the denominator is non-zero, for $r \geq k$, since $\mathbb{P}(R=r)>0$.
A.3. **Verification of (7).** We now show that $\mathbb{E}[\hat{f}|R=r]$ in (5) admits the conditional variance expression (7). We will do so under the following sufficient condition which will permit the interchange of the expectation operator and the infinite sum:

$$ (24) \qquad \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} |\gamma_k||\gamma_l| \frac{1}{\mathbb{P}(R \ge k)} \left| \frac{m}{x_0} - 1 \right|^{k+l} < \infty. $$

We only consider the expression for the second moment since $\mathbb{E}[\mathbb{E}[\hat{f}|R]]^2 = f(m)^2$. Squaring (5) – noting that it is a finite series – and replacing $r$ with the (integer valued) random variable $R$ gives

$$ (25) \qquad \mathbb{E}[\hat{f}|R]^2 = \sum_{k=0}^{\infty} \gamma_k^2 \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)^2} \left(\frac{m}{x_0} - 1\right)^{2k} $$

$$ (26) \qquad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \gamma_k \gamma_l \frac{1\{R \ge l\}}{\mathbb{P}(R \ge k)\mathbb{P}(R \ge l)} \left(\frac{m}{x_0} - 1\right)^{k+l}. $$

Since the terms in the first sum are all positive, Fubini's Theorem for non-negative functions (Billingsley, 1995, Theorem 18.3) allows the interchange of expectation and the sum,

$$ \begin{aligned} \mathbb{E} \left[ \sum_{k=0}^{\infty} \gamma_k^2 \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)^2} \left( \frac{m}{x_0} - 1 \right)^{2k} \right] &= \sum_{k=0}^{\infty} \mathbb{E} \left[ \gamma_k^2 \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)^2} \mathbb{E}[U_{R,k}|R]^2 \right] \\ &= \sum_{k=0}^{\infty} \gamma_k^2 \frac{1}{\mathbb{P}(R \ge k)} \left( \frac{m}{x_0} - 1 \right)^{2k}. \end{aligned} $$

For the second sum, assumption (24) also permits the interchange of expectation and the sum. Thus,

$$ \begin{aligned} \mathbb{E} \left[ \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \gamma_k \gamma_l \frac{1\{R \ge l\}}{\mathbb{P}(R \ge k)\mathbb{P}(R \ge l)} \left( \frac{m}{x_0} - 1 \right)^{k+l} \right] \\ = \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \mathbb{E} \left[ \gamma_k \gamma_l \frac{1\{R \ge l\}}{\mathbb{P}(R \ge k)\mathbb{P}(R \ge l)} \left( \frac{m}{x_0} - 1 \right)^{k+l} \right] \end{aligned} $$

This completes the verification of (7).

Finally, condition (24) holds when the distribution of $R$ satisfies certain conditions. For example, precisely under all the conditions of Proposition 4.
A.4. **Derivation of (9).** We start by computing $\mathbb{E}[(\hat{f})^2|R=r]$:

$$
\begin{align*}
\mathbb{E}[(\hat{f})^2|R=r] &= \mathbb{E}\left[\left(\sum_{k=0}^{\infty} \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)} \gamma_k U_{R,k}\right)^2|R=r\right] \\
&= \mathbb{E}\left[\sum_{k=0}^{\infty} \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)^2} \gamma_k^2 U_{R,k}^2 |R=r\right] \\
&\quad + \mathbb{E}\left[2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{1\{R \ge k\} 1\{R \ge l\}}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \gamma_k \gamma_l U_{R,k} U_{R,l} |R=r\right] \\
&= \sum_{k=0}^{\infty} \frac{1\{R \ge k\}}{\mathbb{P}(R \ge k)^2} \gamma_k^2 U_{R,k}^2 + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{1\{R \ge l\}}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \gamma_k \gamma_l U_{R,k} U_{R,l}.
\end{align*}
$$

This and (25) give

$$
\begin{align*}
\mathrm{var}[\hat{f}|R=r] &= \sum_{k=0}^{\infty} \frac{1\{r \ge k\}}{\mathbb{P}(R \ge k)^2} \gamma_k^2 \left( U_{r,k}^2 - \left(\frac{m}{x_0} - 1\right)^{2k} \right) \\
&\quad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{1\{r \ge l\}}{\mathbb{P}(R \ge k)\mathbb{P}(R \ge l)} \gamma_k \gamma_l \left( U_{r,k} U_{r,l} - \left(\frac{m}{x_0} - 1\right)^{k+l} \right) \\
&= \sum_{k=0}^{\infty} \frac{1\{r \ge k\}}{\mathbb{P}(R \ge k)^2} \gamma_k^2 \mathrm{var}(U_{r,k}) \\
&\quad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{1\{r \ge l\}}{\mathbb{P}(R \ge k)\mathbb{P}(R \ge l)} \gamma_k \gamma_l \mathrm{cov}(U_{r,k}, U_{r,l}).
\end{align*}
$$

Taking expectations w.r.t. R gives the result.

APPENDIX B. PROOFS OF SECTION 2

B.1. **Proof of Proposition 4.** Let $\tilde{f} = \mathbb{E}[\hat{f}|R]$. The lower bound can be found by the law of total variance: $\mathrm{var}[\tilde{f}] \ge \mathrm{var}[\mathbb{E}[\hat{f}|1\{R>0\}]]$ and $\mathbb{E}[\hat{f}|1\{R>0\}] = \gamma_0 + \frac{1\{R>0\}}{\mathbb{P}(R>0)}(f(m) - \gamma_0)$. In the case of $R \sim \text{Geometric}(p)$ we have

$$
\begin{align*}
\mathrm{var}[\tilde{f}] &\ge (f(m) - \gamma_0)^2 \left( \frac{1}{\mathbb{P}(R > 0)} - 1 \right) \\
&= (f(m) - \gamma_0)^2 \left( \frac{1}{1-p} - 1 \right).
\end{align*}
$$

For the upper bound, under Assumption 3, the first sum (7) can be bounded as
follows:

$$
\begin{align*}
\sum_{k=0}^{\infty} \gamma_k^2 \beta_0^{2k} \left( \frac{1}{\mathbb{P}(R \ge k)} - 1 \right) &\le c^2 \sum_{k=0}^{\infty} \left[ \left( \frac{\beta_0^2}{1-p} \right)^k - \beta_0^{2k} \right] \\
&= c^2 \left( \frac{1-p}{1-p-\beta_0^2} - \frac{1}{1-\beta_0^2} \right) \\
&= \frac{c^2 p \beta_0^2}{(1-\beta_0^2)(1-p-\beta_0^2)},
\end{align*}
$$
where the convergence of the infinite sum is guaranteed by the fact that $\beta_0 < 1$ and
$1 - p > \beta_0^2$.

The absolute value of the second term may be bounded in a similar manner:

$$
\begin{align*}
2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} |\gamma_k \gamma_l| \beta_0^{k+l} \left( \frac{1}{\mathbb{P}(R \ge k)} - 1 \right) &\le 2c^2 \sum_{k=0}^{\infty} \left( \frac{1}{\mathbb{P}(R \ge k)} - 1 \right) \sum_{l=k+1}^{\infty} \beta_0^{k+l} \\
&= \frac{2c^2 \beta_0}{1 - \beta_0} \sum_{k=0}^{\infty} \beta_0^{2k} \left( \frac{1}{\mathbb{P}(R \ge k)} - 1 \right) \\
&= \frac{2\beta_0}{1 - \beta_0} \times \frac{c^2 p \beta_0^2}{(1 - \beta_0^2)(1 - p - \beta_0^2)}.
\end{align*}
$$

Combining the two inequalities leads to:

$$
\begin{align*}
\operatorname{var}[\mathbb{E}[\hat{f}|R]] &\le \left(1 + \frac{2\beta_0}{1-\beta_0}\right) \frac{c^2 p \beta_0^2}{(1-\beta_0^2)(1-p-\beta_0^2)} \\
&\le \frac{c^2 \beta_0^2}{(1-\beta_0)^2} \times \frac{p}{1-p-\beta_0^2}.
\end{align*}
$$

B.2. **Proof of Proposition 7.**

B.2.1. *Proof of* $\hat{f}^S \in \mathcal{L}_2$. For brevity, let $\hat{f}$ denote $\hat{f}^S$. By the monotone convergence theorem, $\mathbb{E}[\hat{f}^2] = \sum_{r=0}^{\infty} \mathbb{E}[\hat{f}^2 1\{R=r\}]$. By definition of $\hat{f}$ (see (2)), $\hat{f}1\{R=r\} = \sum_{k=0}^r \gamma_k U_{r,k} 1\{R=r\}/\mathbb{P}(R \ge k)$. Thus,

$$
\begin{align*}
\hat{f}^2 1\{R=r\} &= \sum_{k=0}^{r} \frac{\gamma_k^2 U_{r,k}^2}{\mathbb{P}(R \ge k)^2} 1\{R=r\} + 2 \sum_{k=0}^{r-1} \sum_{l=k+1}^{r} 1\{R=r\} \frac{\gamma_k U_{r,k}}{\mathbb{P}(R \ge k)} \frac{\gamma_l U_{r,l}}{\mathbb{P}(R \ge l)} \\
\mathbb{E}[\hat{f}^2 1\{R=r\}] &= \sum_{k=0}^{r} \frac{\gamma_k^2 \mathbb{E}[U_{r,k}^2]}{\mathbb{P}(R \ge k)^2} \mathbb{E}[1\{R=r\}] + 2 \sum_{k=0}^{r-1} \sum_{l=k+1}^{r} \frac{\gamma_k \gamma_l \mathbb{E}[1\{R=r\}]}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \mathbb{E}[U_{r,k} U_{r,l}].
\end{align*}
$$

The independence of $R$ and $U_{r,k}$ has been invoked in the final expression for
$\mathbb{E}[\hat{f}^2 1\{R=r\}]$. For the simple estimate $\hat{f}^S$,

$$
\mathbb{E}[(U_{r,k})^2] = \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^k = \beta^{2k}.
$$

Along with the Cauchy-Schwarz inequality, gives, $|\mathbb{E}[U_{r,k}U_{r,l}]| \le \mathbb{E}[U_{r,k}^2]^{1/2}\mathbb{E}[U_{r,l}^2]^{1/2} \le \beta^{l+k}$. Thus,

$$
\begin{align*}
\mathbb{E}[\hat{f}^2 1\{R=r\}] &\le \mathbb{P}(R=r) \sum_{k=0}^{r} \frac{\gamma_k^2 \beta^{2k}}{\mathbb{P}(R \ge k)^2} + 2\mathbb{P}(R=r) \sum_{k=0}^{r-1} \sum_{l=k+1}^{r} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k)\mathbb{P}(R \ge l)} \beta^{k+l} \\
\sum_{r=1}^{t} \mathbb{E}[\hat{f}^2 1\{R=r\}] &\le \sum_{r=1}^{t} \mathbb{P}(R=r) \sum_{k=0}^{r} \frac{\gamma_k^2 \beta^{2k}}{\mathbb{P}(R \ge k)^2} \\
&\quad + \sum_{r=1}^{t} 2\mathbb{P}(R=r) \sum_{k=0}^{r-1} \sum_{l=k+1}^{r} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k)\mathbb{P}(R \ge l)} \beta^{k+l}
\end{align*}
$$
The first sum can be expressed as

$$
\begin{align*}
\sum_{r=0}^{t} \mathbb{P}(R=r) \sum_{k=0}^{r} \frac{\gamma_k^2 \beta^{2k}}{\mathbb{P}(R \ge k)^2} &= \sum_{k=0}^{t} \frac{\gamma_k^2 \beta^{2k}}{\mathbb{P}(R \ge k)^2} \sum_{r=k}^{t} \mathbb{P}(R=r) \\
&\le \sum_{k=0}^{t} \frac{\gamma_k^2 \beta^{2k}}{\mathbb{P}(R \ge k)} \\
&\le \sum_{k=0}^{\infty} \frac{c^2 \beta^{2k}}{(1-p)^k} \\
&= c^2 \frac{1-p}{1-p-\beta^2},
\end{align*}
$$

where the final line uses condition $1 - p > \beta^2$. Similarly, the double sum in the right-hand side of $\sum_{r=1}^t \mathbb{E}[\hat{f}^2 1\{R=r\}]$ can be expressed, and then upper bounded as follows:

$$
\begin{align*}
\sum_{r=1}^{t} 2\mathbb{P}(R=r) \sum_{k=0}^{r-1} \sum_{l=k+1}^{r} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \beta^{k+l} \\
&= \sum_{k=0}^{t-1} \sum_{l=k+1}^{t} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \beta^{k+l} \sum_{r=l}^{t} 2\mathbb{P}(R=r) \\
&\le \sum_{k=0}^{t-1} \sum_{l=k+1}^{t} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k)} \beta^{k+l} 2 \\
&\le 2c^2 \sum_{k=0}^{t-1} \sum_{l=k+1}^{t} \frac{\beta^{k+l}}{\mathbb{P}(R \ge k)} \\
&\le 2c^2 \frac{\beta}{1-\beta} \frac{1-p}{1-p-\beta^2}.
\end{align*}
$$

Thus, $\sum_{r=1}^{t} \mathbb{E}[\hat{f}^2 1\{R=r\}]$ is bounded uniformly in $t$ by

$$
(27) \qquad c^2 \frac{1-p}{1-p-\beta^2} + 2c^2 \frac{\beta}{1-\beta} \frac{1-p}{1-p-\beta^2},
$$

and so $\mathbb{E}[\hat{f}^2] < \infty$.

B.2.2. *Proof of (11).* We may bound (9) by the second moment,

$$
\mathbb{E}[\mathrm{var}[\hat{f}|R]] = \mathbb{E}[\mathbb{E}[\hat{f}^2|R]] - \mathbb{E}[\mathbb{E}[\hat{f}|R]^2] \le \mathbb{E}[\hat{f}^2].
$$

Using (27),

$$
\begin{align*}
\mathbb{E}[\hat{f}^2] &\le c^2 \left(1 + \frac{2\beta}{1-\beta}\right) \frac{1-p}{1-p-\beta^2} \\
&= c^2 \frac{1+\beta}{1-\beta} \frac{1-p}{1-p-\beta^2}.
\end{align*}
$$
### B.3. Verification of Example 8.
Recall that by the law of total probability we have

$$ \mathrm{var}[\hat{f}^S] = \mathrm{var}\left[\mathbb{E}[\hat{f}^S|R]\right] + \mathbb{E}\left[\mathrm{var}[\hat{f}^S|R]\right]. $$

If we let $p \to 0$, Proposition 4 guarantees that $\mathrm{var}[\mathbb{E}[\hat{f}^S|R]] \to 0$ regardless of the choice of $f, m, x_0$. The behaviour of $\mathbb{E}[\mathrm{var}[\hat{f}^S|R]]$ as $p$ tends to 0 can be obtained by invoking the dominated convergence theorem to interchange $\lim_{p \to 0}$ and the sums (a result also known as Tannery's Theorem; Boas Jr, 1965):

$$ \begin{aligned} \mathbb{E}[\mathrm{var}[\hat{f}^S|R]] &= \sum_{k=0}^{\infty} \frac{\gamma_k^2}{\mathbb{P}(R \ge k)^2} \left\{ \sum_{r=k}^{\infty} \mathbb{P}(R=r) \mathrm{var}(U_{r,k}^S) \right\} \\ &\quad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \left\{ \sum_{r=l}^{\infty} \mathbb{P}(R=r) \mathrm{cov}(U_{r,k}^S, U_{r,l}^S) \right\} \\ &= \sum_{k=0}^{\infty} \frac{\gamma_k^2}{\mathbb{P}(R \ge k)} \left\{ \left( \frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2 \right)^k - \left(\frac{m}{x_0} - 1\right)^{2k} \right\} \\ &\quad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k)} \left\{ \left( \frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2 \right)^k \left(\frac{m}{x_0} - 1\right)^{l-k} - \left(\frac{m}{x_0} - 1\right)^{k+l} \right\} \\ &\to \sum_{k=0}^{\infty} \gamma_k^2 \left\{ \left( \frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2 \right)^k - \left(\frac{m}{x_0} - 1\right)^{2k} \right\} \\ &\quad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \gamma_k \gamma_l \left\{ \left( \frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2 \right)^k \left(\frac{m}{x_0} - 1\right)^{l-k} - \left(\frac{m}{x_0} - 1\right)^{k+l} \right\}. \end{aligned} $$

For $f(x) = 1/x$, $\gamma_k = (-1)^{k-1}$, and substituting into the above expression gives

$$ \begin{aligned} x_0 \lim_{p \to 0} \mathbb{E}[\mathrm{var}[\hat{f}^S|R]] &= \sum_{k=0}^{\infty} \{\beta^{2k} - \beta_0^{2k}\} \\ &\quad + 2 \sum_{k=0}^{\infty} (-1)^{k-1} \left\{ \beta^{2k} \left(\frac{m}{x_0} - 1\right)^{-k} - \left(\frac{m}{x_0} - 1\right)^k \right\} \sum_{l=k+1}^{\infty} (-1)^{l-1} \left(\frac{m}{x_0} - 1\right)^l \\ &= \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} + 2 \left(\frac{x_0}{m} - 1\right) \sum_{k=0}^{\infty} \left\{ \beta^{2k} - \left(\frac{m}{x_0} - 1\right)^{2k} \right\} \\ &= \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} + 2 \left(\frac{x_0}{m} - 1\right) \left\{ \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} \right\} \\ &= \left(\frac{2x_0}{m} - 1\right) \left\{ \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} \right\}. \end{aligned} $$

### B.4. Proof of Proposition 9.

#### B.4.1. Second moment of cycling variables.
We first show that the assumptions of Proposition 9 imply $\hat{f}^C$ is integrable, and thus $\mathbb{E}[\hat{f}^C|R]$ is well-defined, by appealing to Lemma 2. To use Lemma 2, we need to show Assumption 1 holds for the $U_{r,k}^C$'s.
Let $\tilde{X}_i := (X_i/x_0 - 1)$, then (for $r \ge k$),

$$ \mathbb{E}[|U_{r,k}^{\mathrm{C}}|] \le \mathbb{E}\left[\prod_{i=1}^{k} |\tilde{X}_i|\right] = \left\{\mathbb{E}[|\tilde{X}_i|]\right\}^k \le \left\{\mathbb{E}[|\tilde{X}_i|^2]^{1/2}\right\}^k \le \beta^k, $$

where $\beta$, defined in (10), is the square-root of the second moment of $\tilde{X}_i$. Assumption 3-(i) assumes $|\gamma_k|$ is uniformly bounded by some constant $c$ while condition (10) assumes $\beta < 1$. Hence $\sum |\gamma_k|a_k$ (of Assumption 1) is bounded above by $\sum c\beta_k < \infty$. Thus, Lemma 2 applies to $\hat{f} = \hat{f}^C$.

Before proceeding with the proof of Proposition 9, we derive an expression for the moments of $U_{r,k}^C$. It is convenient to define $\rho := 1 + \sigma^2/(m - x_0)^2$ and observe that $\rho > 1$.

The second moment of the cycling variables can be obtained by direct multiplication, distinguishing between three different cases:

**Lemma 12.** For any $r \ge l \ge k$ we have

(1) If $r \ge l + k$

$$ \mathbb{E}[U_{r,k}^{\mathrm{C}} U_{r,l}^{\mathrm{C}}] = \left(\frac{m}{x_0} - 1\right)^{l+k} \frac{1}{r} \left[ (r - l - k + 1) + 2 \sum_{j=1}^{k-1} \rho^j + (l - k + 1)\rho^k \right], $$

with the convention that $\sum_{j=n_1}^{n_2} = 0$ if $n_1 < n_2$.

(2) If $r < l + k$

$$ \mathbb{E}[U_{r,k}^{\mathrm{C}} U_{r,l}^{\mathrm{C}}] = \left(\frac{m}{x_0} - 1\right)^{l+k} \frac{1}{r} \left[ 2 \sum_{j=l+k-r+1}^{k-1} \rho^j + (l-k+1)\rho^k + (l+k-r+1)\rho^{l+k-r} \right], $$

with the convention that $\sum_{j=n_1}^{n_2} = 0$ if $n_1 < n_2$.

(3) If $r = l$, the cycling variable $U_{r,l}^C$ coincides with $U_{r,l}^S$ and we have

$$ \mathbb{E}[U_{r,k}^{\mathrm{C}} U_{r,l}^{\mathrm{C}}] = \left( \frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2 \right)^k \left(\frac{m}{x_0} - 1\right)^{l-k} = \left(\frac{m}{x_0} - 1\right)^{l+k} \rho^k. $$

*Proof.* We only prove case 1. The proof for the second and third case follows the same lines.

Recall that $\tilde{X}_i = (X_i/x_0 - 1)$, let $\tilde{m} := \mathbb{E}[X_i/x_0 - 1]$ and $\tilde{\sigma}^2 := \mathrm{var}(X_i/x_0 - 1) = \sigma^2/x_0^2$. Expanding $U_{r,l}^C$ we find

$$ \begin{aligned} r \times U_{r,l}^{\mathrm{C}} = & \prod_{i=1}^{l} \tilde{X}_i + \cdots + \prod_{i=k}^{l+k-1} \tilde{X}_i \\ & + \prod_{i=k+1}^{l+k} \tilde{X}_i + \prod_{i=k+2}^{l+k+1} \tilde{X}_i + \cdots + \prod_{i=r-l+1}^{r} \tilde{X}_i \\ & + \left( \prod_{i=r-l+2}^{r} \tilde{X}_i \right) \tilde{X}_1 + \cdots + \tilde{X}_r \prod_{i=1}^{l-1} \tilde{X}_i. \end{aligned} $$

The expansion has been organised in three separate rows for the following reason. Note that each row is non-empty since $r \ge k + l$. (i) We commence by calculating $\mathbb{E}[U_{r,l}^C \prod_{j=1}^k \tilde{X}_j]$. The first and third row involve non-independent products with
the term $\prod_{j=1}^{k} \tilde{X}_j$, whereas for the middle row, the expected value is the product of the expectations. (ii) The previous calculation is sufficient since, by symmetry, $\mathbb{E}[U_{r,k}^{\mathrm{C}} U_{r,l}^{\mathrm{C}}] = \mathbb{E}[U_{r,l}^{\mathrm{C}} \prod_{j=1}^{k} \tilde{X}_j]$.

$$ 
\begin{aligned}
& r \times U_{r,l}^{\mathrm{C}} \times \prod_{j=1}^{k} \tilde{X}_j \\
&= \prod_{j=1}^{k} \tilde{X}_j^2 \prod_{i=k+1}^{l} \tilde{X}_i + \cdots + \tilde{X}_k^2 \prod_{j=1}^{k-1} \tilde{X}_j \prod_{i=k+1}^{l+k-1} \tilde{X}_i \\
&\quad + \prod_{j=1}^{k} \tilde{X}_j \prod_{i=k+1}^{l+k} \tilde{X}_i + \cdots + \prod_{j=1}^{k} \tilde{X}_j \prod_{i=r-l+1}^{r} \tilde{X}_i \\
&\quad + \tilde{X}_1^2 \prod_{j=2}^{k} \tilde{X}_j \prod_{i=r-l+2}^{r} \tilde{X}_i + \cdots + \prod_{j=1}^{k} \tilde{X}_j^2 \prod_{i=r-l+k+1}^{r} \tilde{X}_i + \cdots + \left( \prod_{j=1}^{k} \tilde{X}_j^2 \right) \tilde{X}_r \prod_{i=k+1}^{l-1} \tilde{X}_i,
\end{aligned}
$$

and, taking expectations,

$$ 
\begin{aligned}
\mathbb{E} \left[ r U_{r,l}^{\mathrm{C}} \prod_{j=1}^{k} \tilde{X}_j \right] = & (\tilde{\sigma}^2 + \tilde{m}^2)^k \tilde{m}^{l-k} + \cdots + (\tilde{\sigma}^2 + \tilde{m}^2) \tilde{m}^{l+k-2} \\
& + \tilde{m}^{l+k} + \cdots + \tilde{m}^{l+k} \\
& + (\tilde{\sigma}^2 + \tilde{m}^2) \tilde{m}^{l+k-2} + \cdots + (\tilde{\sigma}^2 + \tilde{m}^2)^k \tilde{m}^{l-k} + \cdots + (\tilde{\sigma}^2 + \tilde{m}^2)^k \tilde{m}^{l-k}.
\end{aligned}
$$

Counting how many times each term in the sum above appears we find that

$$ 
\begin{aligned}
\mathbb{E} \left[ r U_{r,l}^{\mathrm{C}} \prod_{j=1}^{k} \tilde{X}_j \right] &= \tilde{m}^{l+k} (r - (l+k) + 1) + 2 \sum_{j=1}^{k-1} (\tilde{\sigma}^2 + \tilde{m}^2)^j \tilde{m}^{l+k-2j} + (\tilde{\sigma}^2 + \tilde{m}^2)^k \tilde{m}^{l-k} (l-k+1) \\
&= \tilde{m}^{l+k} \left[ (r-l-k+1) + 2 \sum_{j=1}^{k-1} \left( \frac{\tilde{\sigma}^2}{\tilde{m}^2} + 1 \right)^j + (l-k+1) \left( \frac{\tilde{\sigma}^2}{\tilde{m}^2} + 1 \right)^k \right].
\end{aligned}
$$

We can then exploit the symmetry of $U_{r,k}^{\mathrm{C}}$ to conclude

$$ 
\mathbb{E}[U_{r,k}^{\mathrm{C}} U_{r,l}^{\mathrm{C}}] = \frac{\tilde{m}^{l+k}}{r} \left[ (r-l-k+1) + 2 \sum_{j=1}^{k-1} \left( \frac{\tilde{\sigma}^2}{\tilde{m}^2} + 1 \right)^j + (l-k+1) \left( \frac{\tilde{\sigma}^2}{\tilde{m}^2} + 1 \right)^k \right].
$$

Using the fact that $\tilde{m} = (m/x_0 - 1)$ and $\tilde{\sigma}^2 = \sigma^2/x_0^2$ and recalling the definition of $\rho = 1 + \sigma^2/(m - x_0)^2$ we finally obtain

$$ 
\mathbb{E}[U_{r,k}^{\mathrm{C}} U_{r,l}^{\mathrm{C}}] = \left(\frac{m}{x_0} - 1\right)^{l+k} \frac{1}{r} \left[ (r - l - k + 1) + 2 \sum_{j=1}^{k-1} \rho^j + (l - k + 1) \rho^k \right].
$$

□

### B.4.2. Bound on second moments of cycling variables.
To characterise the behaviour of the cycling estimator as the number of terms in the Taylor expansion increases, we establish the following bound on $\mathrm{cov}[U_{r,k}^{\mathrm{C}} U_{r,l}^{\mathrm{C}}]$.
**Lemma 13.** Let $r \ge l \ge k$. Then,

$$ \operatorname{cov}[U_{r,k}^{C}, U_{r,l}^{C}] < \left|\frac{m}{x_0} - 1\right|^{l+k} \times \begin{cases} \rho^k(l+k)/r, & \text{for } r \ge l+k \\ (\rho^k+1), & \text{for } r < l+k \\ (\rho^k-1), & \text{for } r = l. \end{cases} $$

*Proof.* The proofs relies on bounding each one of the second moments obtained in Lemma 12 separately. Recall that $\operatorname{cov}[U_{r,k}^C U_{r,l}^C] = \mathbb{E}[U_{r,k}^C U_{r,l}^C] - \mathbb{E}[U_{r,k}^C] \mathbb{E}[U_{r,l}^C]$. Using the expression for $\mathbb{E}[U_{r,k}^C U_{r,l}^C]$ in Lemma 12-1, since $\rho > 1$, we have

$$ \begin{aligned} \operatorname{cov}[U_{r,k}^C U_{r,l}^C] &= \left(\frac{m}{x_0} - 1\right)^{l+k} \frac{1}{r} \left[ (r - l - k + 1) + 2 \sum_{j=1}^{k-1} \rho^j + (l - k + 1)\rho^k - r \right] \\ &< \left|\frac{m}{x_0} - 1\right|^{l+k} \frac{1}{r} \left[ 1 + 2 \sum_{j=1}^{k-1} \rho^j + (l - k + 1)\rho^k \right] \\ &< \frac{\rho^k}{r} \left|\frac{m}{x_0} - 1\right|^{l+k} \left[ 1 + 2 \sum_{j=1}^{k-1} \rho^{j-k} + (l - k + 1) \right] \\ &< \frac{\rho^k}{r} \left|\frac{m}{x_0} - 1\right|^{l+k} \left[ 1 + 2(k - 1) + (l - k + 1) \right] \\ &= \frac{\rho^k}{r} (l + k) \left|\frac{m}{x_0} - 1\right|^{l+k}. \end{aligned} $$

Similarly, for case 2 we find

$$ \begin{aligned} \mathbb{E}[U_{r,k}^C U_{r,l}^C] &= \frac{1}{r} \left(\frac{m}{x_0} - 1\right)^{l+k} \rho^k \left[ 2 \sum_{j=l+k-r+1}^{k-1} \rho^{j-k} + (l-k+1) + (l+k-r+1)\rho^{l-r} \right] \\ &< \frac{1}{r} \left|\frac{m}{x_0} - 1\right|^{l+k} \rho^k \left[ 2(r-l-1) + (l-k+1) + (l+k-r+1) \right] \\ &= \left|\frac{m}{x_0} - 1\right|^{l+k} \rho^k, \end{aligned} $$

where we used the fact that $\rho > 1$. It follows that

$$ \operatorname{cov}[U_{r,k}^C U_{r,l}^C] < \left|\frac{m}{x_0} - 1\right|^{l+k} \left( \rho^k - \operatorname{sgn}\left(\left(\frac{m}{x_0} - 1\right)^{l+k}\right) \right) \le \left|\frac{m}{x_0} - 1\right|^{l+k} (\rho^k + 1). $$

Finally, case 3 follows straightforwardly

$$ \operatorname{cov}[U_{r,k}^C U_{r,l}^C] = \left(\frac{m}{x_0} - 1\right)^{l+k} (\rho^k - 1). $$

Collecting all the stated, for cases 1, 2 and 3, completes the proof. □

B.4.3. *Proof of Proposition 9.* Before proceeding to the proof of Proposition 9 we state and prove the following simple result:
**Lemma 14 (Expectation of Reciprocal).** If $R \sim \text{Geometric}(p)$ then for all $k \ge 1$,

$$ \mathbb{E} \left[ \frac{1}{R} \mid R \ge k \right] = \frac{1}{\mathbb{P}(R \ge k)} \sum_{r=k}^{\infty} \frac{1}{r} \mathbb{P}(R=r) \le \frac{p \log(1/p)}{1-p}. $$

*Proof.* By definition we have

$$ \begin{aligned} \frac{1}{\mathbb{P}(R \ge k)} \sum_{r=k}^{\infty} \frac{1}{r} \mathbb{P}(R=r) &= p \sum_{r=k}^{\infty} \frac{(1-p)^{r-k}}{r} \\ &= p \sum_{h=1}^{\infty} \frac{(1-p)^{h-1}}{h+k-1} \\ &\le \frac{p}{1-p} \sum_{h=1}^{\infty} \frac{(1-p)^h}{h} \\ &= \frac{p}{1-p} \log(1/p). \end{aligned} $$

□

We are now ready to prove Proposition 9. First, we observe that for $k = 0$ we have $U_{r,0}^C = 1$, thus $\text{var}(U_{r,0}^C) = 0$ and we have

$$ \mathbb{E} [\text{var}[\hat{f}^C|R]] = \sum_{k=1}^{\infty} \frac{\gamma_k^2}{\mathbb{P}(R \ge k)^2} \left\{ \sum_{r=k}^{\infty} \mathbb{P}(R=r) \text{var}(U_{r,k}^C) \right\} \\ + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \left\{ \sum_{r=l}^{\infty} \mathbb{P}(R=r) \text{cov}(U_{r,k}^C, U_{r,l}^C) \right\}. $$

We use Lemma 13 to bound the terms $\text{cov}(U_{r,k}^C, U_{r,l}^C)$. For the case $r < l + k$, we use the simplified bound $\rho^k + 1 \le 2\rho^k$ and then $\rho^k < \rho^k(l+k)/r \le \rho^k(2l)/r$. When $r = l$, clearly $\rho^k - 1 < \rho^k(2l)/r$.

Combining the above bounds from Lemma 13 with Lemma 14, gives

$$ \begin{aligned} \mathbb{E} [\text{var}[\hat{f}^C|R]] &\le \sum_{k=1}^{\infty} \frac{\gamma_k^2}{\mathbb{P}(R \ge k)^2} \left\{ 2k \left| \frac{m}{x_0} - 1 \right|^{2k} \rho^k \sum_{r=k}^{\infty} \frac{1}{r} \mathbb{P}(R=r) \right\} \\ &\quad + 2 \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{|\gamma_k \gamma_l|}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \left\{ 4l \left| \frac{m}{x_0} - 1 \right|^{l+k} \rho^k \sum_{r=l}^{\infty} \frac{1}{r} \mathbb{P}(R=r) \right\} \\ &\le \frac{p \log(1/p)}{1-p} \sum_{k=1}^{\infty} \frac{\gamma_k^2}{\mathbb{P}(R \ge k)} 2k \left| \frac{m}{x_0} - 1 \right|^{2k} \rho^k \\ &\quad + 2 \frac{p \log(1/p)}{1-p} \sum_{k=0}^{\infty} \sum_{l=k+1}^{\infty} \frac{|\gamma_k \gamma_l|}{\mathbb{P}(R \ge k)} 4l \left| \frac{m}{x_0} - 1 \right|^{l+k} \rho^k. \end{aligned} $$

Using Assumption 3 and recalling that

$$ \beta^{2k} = \left( \frac{\sigma^2}{x_0^2} + \left( \frac{m}{x_0} - 1 \right)^2 \right)^k = \left( \frac{m}{x_0} - 1 \right)^{2k} \rho^k, $$
that $\beta_0 = |m/x_0 - 1| < 1$, we have

$$
\begin{align*}
\mathbb{E}[\mathrm{var}[\hat{f}^C|R]] &\le 2c^2 \frac{p \log(1/p)}{1-p} \left\{ \sum_{k=1}^{\infty} k \frac{\beta^{2k}}{(1-p)^k} + 4 \sum_{k=0}^{\infty} \frac{\beta_0^k \rho^k}{(1-p)^k} \sum_{l=k+1}^{\infty} l \beta_0^l \right\} \\
&= \frac{2c^2 p \log(1/p)}{1-p} \left\{ \frac{\beta^2(1-p)}{(1-p-\beta^2)^2} + 4 \frac{\beta_0}{(1-\beta_0)^2} \sum_{k=0}^{\infty} \frac{\beta_0^{2k} \rho^k}{(1-p)^k} (k(1-\beta_0) + 1) \right\} \\
&\le \frac{2c^2 p \log(1/p)}{1-p} \left\{ \frac{\beta^2(1-p)}{(1-p-\beta^2)^2} + 4 \frac{\beta_0}{(1-\beta_0)^2} \sum_{k=0}^{\infty} \frac{\beta^{2k}}{(1-p)^k} (k+1) \right\} \\
&= 2c^2 \frac{p \log(1/p)}{(1-p-\beta^2)^2} \left( \beta^2 + 4 \frac{\beta_0(1-p)}{(1-\beta_0)^2} \right),
\end{align*}
$$

where we used that $\sum_{k=1}^{\infty} kg^k = g/(1-g)^2$ and $\sum_{k=0}^{\infty} (k+1)g^k = 1/(1-g)^2$ for $g \in (0,1)$.

B.5. **Proof of Proposition 11.** The proof of the first assertion, that $\widehat{\nabla_{\theta} f}^S$ and $\widehat{\nabla_{\theta} f}^C$ have finite second moments, follows closely the steps of Proposition 7's proof (that $\hat{f}^S$ has a finite second moment.) This part is thus omitted.

For the remaining three statements on the variance, the following variance decomposition of $\widehat{\nabla_{\theta} f}$ will be used,

$$
\mathrm{var}[\widehat{\nabla}_{\theta} f] = \mathrm{var}[\mathbb{E}[\widehat{\nabla}_{\theta} f | R]] + \mathbb{E}[\mathrm{var}[\widehat{\nabla}_{\theta} f | R]] ,
$$

where the first term can be expressed as

$$
(28) \quad \mathrm{var} \left[ \mathbb{E}[\widehat{\nabla}_{\theta} f | R] \right] = \sum_{k=1}^{\infty} k^2 \gamma_k^2 \left( \frac{m}{x_0} - 1 \right)^{2k-2} \frac{\nabla_{\theta} m(\theta)^2}{x_0^2} \left( \frac{1}{\mathbb{P}(R \ge k)} - 1 \right) \\
\qquad + 2 \sum_{k=1}^{\infty} \sum_{l=k+1}^{\infty} kl \gamma_k \gamma_l \left( \frac{m}{x_0} - 1 \right)^{k+l-2} \frac{\nabla_{\theta} m(\theta)^2}{x_0^2} \left( \frac{1}{\mathbb{P}(R \ge k)} - 1 \right).
$$

As for the second term, it is

$$
(29) \quad \mathbb{E}[\mathrm{var}[\widehat{\nabla}_{\theta} f | R]] = \frac{1}{x_0^2} \sum_{k=1}^{\infty} k^2 \frac{\gamma_k^2}{\mathbb{P}(R \ge k)^2} \left\{ \sum_{r=k}^{\infty} \mathbb{P}(R=r) \mathrm{var}(W_{r,k}) \right\} \\
\qquad + \frac{2}{x_0^2} \sum_{k=1}^{\infty} \sum_{l=k+1}^{\infty} kl \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \left\{ \sum_{r=l}^{\infty} \mathbb{P}(R=r) \mathrm{cov}(W_{r,k}, W_{r,l}) \right\}
$$

where the expression of $W_{r,k}$ depends on the chosen estimator, either simple or
cycling, for $(m(\theta)/x_0 - 1)^{k-1} \nabla_{\theta} m(\theta)$.

Statements (1)–(3) of Proposition 11 are established as follows:

(1) The stated result for $\mathrm{var}[\mathbb{E}[\widehat{\nabla}_{\theta} f | R]]$ is established by Proposition 15. Note that $\mathbb{E}[\widehat{\nabla}_{\theta} f | R]$ has the same expression for both the simple and cycling case.

(2) Example 16 shows $\lim_{p \to 0} \mathbb{E}[\mathrm{var}[\widehat{\nabla}_{\theta} f^S | R]] > 0$ for $f(\cdot) = \log(\cdot)$. The necessary supporting results to enable this calculation are presented just before this example.

(3) Proposition 17 verifies the final assertion for the cycle estimate.
**Proposition 15.** Let Assumption 3 hold with $m$ therein replaced by $m(\theta)$. If there exists $p \in (0, 1 - \beta_0^2)$ such that $\mathbb{P}(R \ge k) \ge (1-p)^{k-1}$, then a bound for (28) is

$$ \operatorname{var}[\mathbb{E}[\widehat{\nabla_{\theta} f}|R]] \le c^2 \frac{\nabla_{\theta} m(\theta)^2}{x_0^2} \frac{1}{(1 - \beta_0)^2} \times \left( \frac{(\beta_0^2 + \beta_0^4)(1 - \beta_0)^2 + 4\beta_0}{(1 - \beta_0^2/(1 - p))^3} - \frac{(\beta_0^2/(1 - p) + \beta_0^4/(1 - p)^2)(1 - \beta_0)^2 + 4\beta_0}{(1 - \beta_0^2)^3} \right). $$

It is can be shown that the right-hand side is $O(p)$ as $p \to 0$. The proof of this result is very similar to that of Proposition 4 and is therefore omitted.

B.5.1. *Simple gradient estimator.* Using the definition of $W_{r,k}^S$ in (16) it is easy to see that for any $0 \le k \le l$,

$$ (30) \qquad \operatorname{var}(W_{r,k}^S) = \left( \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 \right)^{k-1} s^2(\theta) - \left( \frac{m(\theta)}{x_0} - 1 \right)^{2k-2} \nabla_{\theta} m(\theta)^2, \\ (31) \qquad \operatorname{cov}(W_{r,k}^S, W_{r,l}^S) = \left( \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 \right)^{k-1} \left( \frac{m(\theta)}{x_0} - 1 \right)^{l-k-1} \nabla_{\theta} m(\theta) t(\theta) $$

$$ \operatorname{cov}(W_{r,k}^S, W_{r,l}^S) = \left( \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 \right)^{k-1} \left( \frac{m(\theta)}{x_0} - 1 \right)^{l-k-1} \nabla_{\theta} m(\theta) t(\theta) - \left( \frac{m(\theta)}{x_0} - 1 \right)^{k+l-2} \nabla_{\theta} m(\theta)^2, $$

where we denoted $s^2(\theta) := \mathbb{E}[G_k^2]$ and $t(\theta) := \mathbb{E}[(X_k/x_0 - 1)G_k]$.

To obtain the limit $\lim_{p \to 0} \mathbb{E}[\operatorname{var}(\widehat{\nabla_{\theta} f}^S | R)]$, substitute (30) and (31) into (29) and use Tannery's Theorem; Boas Jr, 1965 to interchange the limit and the sum (as previously shown for the estimator of $f(m)$):

$$ x_0^2 \mathbb{E} \left[ \operatorname{var}[\widehat{\nabla_{\theta} f}^S | R] \right] \rightarrow \sum_{k=1}^{\infty} k^2 \gamma_k^2 \left\{ \left( \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 \right)^{k-1} s^2(\theta) - \left( \frac{m(\theta)}{x_0} - 1 \right)^{2k-2} \nabla_{\theta} m(\theta)^2 \right\} \\ + 2 \sum_{k=1}^{\infty} \sum_{l=k+1}^{\infty} kl\gamma_k\gamma_l \left( \frac{m(\theta)}{x_0} - 1 \right)^{l-k-1} \nabla_{\theta} m(\theta) \\ \times \left\{ \left( \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 \right)^{k-1} t(\theta) - \left( \frac{m(\theta)}{x_0} - 1 \right)^{2k-1} \nabla_{\theta} m(\theta) \right\} $$

The following example shows $\lim_{p \to 0} \mathbb{E}[\operatorname{var}[\widehat{\nabla_{\theta} f}^S | R]]$ can be non-zero.
**Example 16.** Consider $f(x) = \log x$ so that $\gamma_k = (-1)^{k-1}/k$ and assume that we want to estimate $f(m)$ for $m > 0$. Then, we have

$$
\begin{align*}
x_0^2 \lim_{p \to 0} \mathbb{E} [\mathrm{var}[\widehat{\nabla_{\theta} f}^S | R]] &= \sum_{k=1}^{\infty} \{ \beta^{2(k-1)} s^2(\theta) - \beta_0^{2(k-1)} \nabla_{\theta} m(\theta)^2 \} \\
&\quad + 2 \nabla_{\theta} m(\theta) \sum_{k=1}^{\infty} \{ \beta^{2(k-1)} t(\theta) - \beta_0^{2(k-1)} \beta_0 \nabla_{\theta} m(\theta) \} \\
&\quad \times \sum_{l=k+1}^{\infty} (-1)^{k+l-2} \left( \frac{m(\theta)}{x_0} - 1 \right)^{l-k-1} \\
&= \frac{s^2(\theta)}{1 - \beta^2} - \frac{\nabla_{\theta} m(\theta)^2}{1 - \beta_0^2} \\
&\quad - \frac{2 \nabla_{\theta} m(\theta)}{1 + \beta_0} \sum_{k=1}^{\infty} \{ \beta^{2(k-1)} t(\theta) - \beta_0^{2(k-1)} \beta_0 \nabla_{\theta} m(\theta) \} \\
&= \frac{s^2(\theta)}{1 - \beta^2} - \frac{\nabla_{\theta} m(\theta)^2}{1 - \beta_0^2} \\
&\quad - \frac{2 \nabla_{\theta} m(\theta)}{1 + \beta_0} \left( \frac{t(\theta)}{1 - \beta^2} - \frac{\beta_0 \nabla_{\theta} m(\theta)}{1 - \beta_0^2} \right).
\end{align*}
$$

Since $s^2(\theta) = \mathbb{E}[G_k^2] \ge \mathbb{E}[G_k]^2 = \nabla_{\theta} m(\theta)^2$ we can lower bound the first term with

$$
\frac{s^2(\theta)}{1 - \beta^2} - \frac{\nabla_{\theta} m(\theta)^2}{1 - \beta_0^2} \geq \nabla_{\theta} m(\theta)^2 \left( \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} \right).
$$

Then, assume that we can sample $G_k$ and $X_k$ independently, so that $t(\theta) = \beta_0 \nabla_{\theta} m(\theta)$.
It follows that the second term becomes

$$
- \frac{2 \nabla_{\theta} m(\theta)^2 \beta_0}{1 + \beta_0} \left( \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} \right) = 2 \nabla_{\theta} m(\theta)^2 \left( \frac{x_0}{m} - 1 \right) \left( \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} \right).
$$

Hence, we can lower bound

$$
x_0^2 \lim_{p \to 0} \mathbb{E} [\mathrm{var}[\widehat{\nabla_{\theta} f}^S | R]] \geq \nabla_{\theta} m(\theta)^2 \left( \frac{1}{1 - \beta^2} - \frac{1}{1 - \beta_0^2} \right) \left( \frac{2x_0}{m} - 1 \right).
$$

If $\sigma^2 > 0$, we have $\beta^2 > \beta_0^2$ and, since we need to select $x_0 > m/2$ to guarantee
that the Taylor sum (14) converges, we have $\lim_{p \to 0} \mathbb{E}[\mathrm{var}[\widehat{\nabla_{\theta} f}^S | R]] > 0$, showing
that the variance of the simple estimator does not converge to 0 in general.

B.5.2. *Cycling gradient estimator.* For the cycling estimator, under Assumption 3 (with $m$ replaced by $m(\theta)$ therein), $\mathbb{E}[\mathrm{var}(\widehat{\nabla_{\theta} f^C} | R)]$ satisfies a result similar to that in Proposition 9. The proof of this result can established by following the steps outlined to prove Proposition 9 (that gave the corresponding bound for $\mathbb{E}[\mathrm{var}(\hat{f}^C | R)]$.) In particular, substitute into (29) the bounds on the second moments of $W_{r,k}^C$ given in Lemma 19, and then also use Lemma 14.

**Proposition 17.** Let Assumption 3 hold with $m$ therein replaced by $m(\theta)$. Assume further that

$$
\beta^2 := \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 < 1,
$$
and $R \sim \text{Geometric}(p)$ with $p \in (0, 1 - \beta^2)$. Then,

$$ \mathbb{E} \left[ \text{var}(\widehat{\nabla_{\theta} f}^{C} | R) \right] \le \frac{8c^2 s(\theta)^2 p \log(1/p)}{x_0^2 (1-p)} \times \left\{ \frac{\beta^2(1-p)(\beta^2 + 2(1-p))}{\beta_0(1-p-\beta^2)^4} + \frac{\beta_0^2 + 2}{(1-\beta_0)^3} \frac{\beta^2(\beta^4 + 4\beta^2(1-p) + (1-p)^2)}{(1-p-\beta^2)^4} \right\}, $$

*Proof.* Recall that from the law of total variance decomposition we have

$$ (32) \quad \mathbb{E} \left[ \text{var}[\widehat{\nabla_{\theta} f}^{C} | R] \right] = \frac{1}{x_0^2} \sum_{k=1}^{\infty} k^2 \frac{\gamma_k^2}{\mathbb{P}(R \ge k)^2} \left\{ \sum_{r=k}^{\infty} \mathbb{P}(R=r) \text{var}(W_{r,k}^{C}) \right\} \\ + \frac{2}{x_0^2} \sum_{k=1}^{\infty} \sum_{l=k+1}^{\infty} kl \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \left\{ \sum_{r=l}^{\infty} \mathbb{P}(R=r) \text{cov}(W_{r,k}^{C}, W_{r,l}^{C}) \right\}. $$

We use Lemma 19 to bound the terms $\text{cov}(W_{r,k}^C, W_{r,l}^C)$, and Lemma 14 to bound the expectation w.r.t. $R$. Since $\beta_0^2 < \beta^2 < 1$, $\max(\beta_0, 1) = 1$.

$$ \begin{aligned} \mathbb{E} \left[ \text{var}[\widehat{\nabla_{\theta} f}^{C} | R] \right] & \le \frac{2s(\theta)^2}{x_0^2} \sum_{k=1}^{\infty} (2k-2)k^2 \frac{\gamma_k^2}{\mathbb{P}(R \ge k)^2} \rho(\theta)^{k-1} \beta_0^{2k-3} \left\{ \sum_{r=k}^{\infty} \mathbb{P}(R=r) \frac{1}{r} \right\} \\ & \quad + \frac{4s(\theta)^2}{x_0^2} \sum_{k=1}^{\infty} \sum_{l=k+1}^{\infty} (2l-2)kl \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k) \mathbb{P}(R \ge l)} \rho(\theta)^{k-1} \beta_0^{k+l-3} \left\{ \sum_{r=l}^{\infty} \mathbb{P}(R=r) \frac{1}{r} \right\} \\ & \le \frac{p \log(1/p)}{1-p} \frac{2s(\theta)^2}{x_0^2} \sum_{k=1}^{\infty} (2k-2)k^2 \frac{\gamma_k^2}{\mathbb{P}(R \ge k)} \rho(\theta)^{k-1} \beta_0^{2k-3} \\ & \quad + \frac{p \log(1/p)}{1-p} \frac{4s(\theta)^2}{x_0^2} \sum_{k=1}^{\infty} \sum_{l=k+1}^{\infty} (2l-2)kl \frac{\gamma_k \gamma_l}{\mathbb{P}(R \ge k)} \rho(\theta)^{k-1} \beta_0^{k+l-3}. \end{aligned} $$

Using Assumption 3 and recalling that

$$ \beta^{2k} = \left( \frac{\sigma^2(\theta)}{x_0^2} + \left( \frac{m(\theta)}{x_0} - 1 \right)^2 \right)^k = \beta_0^{2k} \rho(\theta)^k, $$
we have

$$
\begin{align*}
& \mathbb{E}[\widehat{\mathrm{var}[\nabla_{\theta} f^C | R]}}] \\
& \le \frac{p \log(1/p)}{1-p} \frac{2c^2 s(\theta)^2}{x_0^2} \sum_{k=1}^{\infty} \frac{(2k-2)k^2}{(1-p)^k} \rho(\theta)^{k-1} \beta_0^{2k-3} \\
& \quad + \frac{p \log(1/p)}{1-p} \frac{4c^2 s(\theta)^2}{x_0^2} \sum_{k=1}^{\infty} \sum_{l=k+1}^{\infty} \frac{(2l-2)kl}{(1-p)^k} \rho(\theta)^{k-1} \beta_0^{k+l-3} \\
& = \frac{p \log(1/p)}{1-p} \frac{2c^2 s(\theta)^2}{\beta_0 x_0^2} \sum_{k=1}^{\infty} (2k-2)k^2 \frac{\beta^{2(k-1)}}{(1-p)^k} \\
& \quad + \frac{p \log(1/p)}{1-p} \frac{4c^2 s(\theta)^2}{x_0^2} \sum_{k=1}^{\infty} \frac{k}{(1-p)^k} \beta_0^{k-3} \rho(\theta)^{k-1} \sum_{l=k+1}^{\infty} (2l-2)l\beta_0^l \\
& = \frac{2c^2 s(\theta)^2 p \log(1/p)}{\beta_0 x_0^2 (1-p)} \frac{4\beta^2(1-p)(\beta^2+2(1-p))}{(1-p-\beta^2)^4} \\
& \quad + \frac{4c^2 s(\theta)^2 p \log(1/p)}{x_0^2 (1-p)} \sum_{k=1}^{\infty} \frac{k}{(1-p)^k} \beta_0^{2k-2} \rho(\theta)^{k-1} \frac{2(\beta_0^2 k^2 - 2k^2 \beta_0 + k^2 - k\beta_0^2 + k + 2\beta_0)}{(1-\beta_0)^3} \\
& \le \frac{2c^2 s(\theta)^2 p \log(1/p)}{\beta_0 x_0^2 (1-p)} \frac{4\beta^2(1-p)(\beta^2+2(1-p))}{(1-p-\beta^2)^4} \\
& \quad + \frac{8c^2 s(\theta)^2 p \log(1/p)}{x_0^2 (1-p)} \frac{\beta_0^2 + 2}{(1-\beta_0)^3} \sum_{k=1}^{\infty} k^3 \frac{\beta^{2(k-1)}}{(1-p)^k} \\
& \le \frac{8c^2 s(\theta)^2 p \log(1/p)}{x_0^2 (1-p)} \\
& \quad \times \left\{ \frac{\beta^2(1-p)(\beta^2+2(1-p))}{\beta_0(1-p-\beta^2)^4} + \frac{\beta_0^2 + 2}{(1-\beta_0)^3} \frac{\beta^2(\beta^4 + 4\beta^2(1-p) + (1-p)^2)}{(1-p-\beta^2)^4} \right\},
\end{align*}
$$

where we used that $\sum_{k=0}^{\infty} 2k(k+1)^2 g^k = 4g(g+2)/(1-g)^4$ and $\sum_{k=0}^{\infty} k^3 g^k = g(g^2+4g+1)/(1-g)^4$. $\square$

B.5.3. *Second moments of cycling variables.* The second moments of the cycling variables $W_{r,k}^C$ can be obtained by a simple counting argument identical to that used in the proof of Lemma 12 as shown in the Proposition below.

**Lemma 18.** Let $s(\theta)^2 = \mathbb{E}[G_1^2]$, $t(\theta) = \mathbb{E}((X_1/x_0 - 1)G_1)$ and

$$
\rho(\theta) = \frac{\mathbb{E}((X_1/x_0 - 1)^2)}{(\mathbb{E}(X_1/x_0 - 1))^2}.
$$
For any $r \ge l \ge k$,

$$
\begin{align}
& \mathbb{E} (W_{r,k}^{\mathrm{C}} W_{r,l}^{\mathrm{C}}) \tag{33} \\
&= \frac{1}{r} \sum_{i=1}^{1} \rho(\theta)^{k-i} \left( \frac{m(\theta)}{x_0} - 1 \right)^{k+l-3} \left( 1\{l > k\}t(\theta)\nabla m(\theta) + 1\{l = k\}s(\theta)^2 \left( \frac{m(\theta)}{x_0} - 1 \right) \right) \nonumber \\
& \qquad + \frac{1}{r} \sum_{i=2, i<k}^{r-l+1} \rho(\theta)^{k-i} \left( \frac{m(\theta)}{x_0} - 1 \right)^{k+l-3} t(\theta)\nabla m(\theta) \tag{34} \\
& \qquad + \frac{1}{r} \sum_{i=1, i \ge k}^{r-l+1} \left( \frac{m(\theta)}{x_0} - 1 \right)^{k+l-2} (\nabla m(\theta))^2 \tag{35} \\
& \qquad + \frac{1}{r} \sum_{i=1}^{k-1} \rho(\theta)^{i-1} \left( \frac{m(\theta)}{x_0} - 1 \right)^{k+l-3} t(\theta)\nabla m(\theta) \tag{36} \\
& \qquad + \frac{1}{r} \sum_{i=k}^{k} 1\{l > k\} \rho(\theta)^{k-1} s(\theta)^2 \left( \frac{m(\theta)}{x_0} - 1 \right)^{l+k-2} \nonumber \\
& \qquad + \frac{1}{r} \sum_{i=k+1}^{l-1} \rho(\theta)^{k-1} \left( \frac{m(\theta)}{x_0} - 1 \right)^{l+k-3} t(\theta)\nabla m(\theta), \tag{38}
\end{align}
$$

where empty sums, $\sum_{i=i_1}^{i_2}$ with $i_2 < i_1$, are by convention zero.

*Proof.* Equivalently, we can also write $\rho(\theta) = 1 + \sigma(\theta)^2/(m(\theta) - x_0)^2$. The sums with single terms, namely (33) and (37), single the unique cases that involve $s(\theta)^2$, as explained more below. For integers $i, j$, with $i < j$, let $i : j$ denote the sequence $i, i+1, \dots, j$. For $r \ge l$, let

$$
(39) \quad W_{r,l}^{\mathrm{C}} = \frac{1}{r} \sum_{i=1}^{r-l+1} V(i: i+l-1) + \frac{1}{r} \sum_{i=1}^{l-1} \bar{V}(r-(l-i)+1:r)V(1:i)
$$

where

$$
V(i:j) = \prod_{k=i}^{j-1} \left( \frac{X_k}{x_0} - 1 \right) \times G_j, \quad j > i
$$

and for $i = j$,

$V(j:j) = G_j.$

Finally,

$$
\bar{V}(i:j) = \prod_{k=i}^{j} \left( \frac{X_k}{x_0} - 1 \right), \quad j \ge i.
$$

For $r \ge l \ge k$, by symmetry of the circular shift that defines $W_{r,l}^C$, we have

$$
(40) \qquad \mathbb{E} (W_{r,k}^{\mathrm{C}} W_{r,l}^{\mathrm{C}}) = \mathbb{E} (V(1:k) \times W_{r,l}^{\mathrm{C}}),
$$
where, using (39), we have
(41)

$$ \frac{1}{r} \sum_{i=1}^{r-l+1} \mathbb{E}(V(1:k)V(i:i+l-1)) + \frac{1}{r} \sum_{i=1}^{l-1} \mathbb{E}(V(1:k)\bar{V}(r-(l-i)+1:r)V(1:i)), $$

Evaluating the products in (41) due to the first sum: For $i > k$

$$ \begin{align} \mathbb{E}(V(1:k)V(i:i+l-1)) &= \mathbb{E}(V(1:k))\mathbb{E}(V(i:i+l-1)) \nonumber \\ &= \left(\frac{m}{x_0} - 1\right)^{k-1} \nabla m \times \left(\frac{m}{x_0} - 1\right)^{l-1} \nabla m \tag{42} \\ &= \left(\frac{m}{x_0} - 1\right)^{k+l-2} (\nabla m)^2, \nonumber \end{align} $$

where the $\theta$ dependence has been omitted for brevity. This accounts for (35). For $1 < i \le k$,

$$ \begin{align*} \mathbb{E}(V(1:k)V(i:i+l-1)) &= \mathbb{E}(V(1:k)\bar{V}(i:k)V((k+1):(i+l-1))) \\ &= \mathbb{E}\left(V(1:k)\left(\frac{X_i}{x_0} - 1\right) \cdots \left(\frac{X_k}{x_0} - 1\right)\right) \mathbb{E}\left(V((k+1):(i+l-1))\right) \\ &= \mathbb{E}\left(V(1:k)\left(\frac{X_i}{x_0} - 1\right) \cdots \left(\frac{X_k}{x_0} - 1\right)\right) \left(\frac{m}{x_0} - 1\right)^{i+l-k-2} \nabla m \end{align*} $$

and

$$ \begin{align*} \mathbb{E}\left(V(1:k)\left(\frac{X_i}{x_0} - 1\right) \cdots \left(\frac{X_k}{x_0} - 1\right)\right) \\ &= \mathbb{E}\left(\left(\frac{X_1}{x_0} - 1\right) \cdots \left(\frac{X_{i-1}}{x_0} - 1\right) \left(\frac{X_i}{x_0} - 1\right)^2 \cdots \left(\frac{X_{k-1}}{x_0} - 1\right)^2 G_k\left(\frac{X_k}{x_0} - 1\right)\right) \\ &= \left(\frac{m}{x_0} - 1\right)^{i-1} \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{k-i} t. \end{align*} $$

Thus

$$ \begin{align} \mathbb{E}(V(1:k)V(i:i+l-1)), \quad 1 < i \le k, \tag{43} \\ &= \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{k-i} \left(\frac{m}{x_0} - 1\right)^{l-k+2i-3} \times t \times \nabla m \nonumber \\ &= \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{k-i} \left(\frac{m}{x_0} - 1\right)^{-2(k-i)} \left(\frac{m}{x_0} - 1\right)^{l+k-3} t \nabla m \nonumber \\ &= \rho^{k-i} \left(\frac{m}{x_0} - 1\right)^{k+l-3} t \nabla m, \nonumber \end{align} $$

which accounts for (34).
For $i = 1$, if $k = l$,

$$ (44) \qquad \begin{aligned} \mathbb{E}(V(1:k)V(1:l)) &= \mathbb{E}(V(1:k)^2) \\ &= \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{k-1} s^2 \\ &= \rho^{k-1} \left(\frac{m}{x_0} - 1\right)^{2(k-1)} s^2, \end{aligned} $$

which accounts for (33) when $\mathbb{1}\{l=k\} = 1$. If $l > k$,

$$ (45) \qquad \begin{aligned} \mathbb{E}(V(1:k)V(1:l)) &= \mathbb{E}(V(1:k)\bar{V}(1:k)V((k+1):l)) \\ &= \rho^{k-i} \left(\frac{m}{x_0} - 1\right)^{k+l-3} t\nabla m. \end{aligned} $$

Note this expected value can be studied as in the case $1 < i \le k$ above to yield the same answer as in (43). Thus accounting for (33) when $\mathbb{1}\{l > k\} = 1$.

Evaluating the products in (41) due to the second sum of $W_{r,l}^C$ in (39): For $1 \le i < k$,

$$ V(1:k)\bar{V}(r-(l-i)+1:r)V(1:i) = \bar{V}(1:i)V(i+1:k)\bar{V}(r-(l-i)+1:r)V(1:i) $$

and thus,

$$ (46) \qquad \begin{aligned} \mathbb{E}(\bar{V}(1:i)V(i+1:k)\bar{V}(r-(l-i)+1:r)V(1:i)) \\ &= \mathbb{E}(\bar{V}(1:i)V(1:i))\mathbb{E}(V(i+1:k))\mathbb{E}(\bar{V}(r-(l-i)+1:r)) \\ &= \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{i-1} t \left(\frac{m}{x_0} - 1\right)^{k-i-1} \nabla m \left(\frac{m}{x_0} - 1\right)^{l-i} \\ &= \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{i-1} \left(\frac{m}{x_0} - 1\right)^{k+l-2i-1} t \nabla m \\ &= \rho^{i-1} \left(\frac{m}{x_0} - 1\right)^{k+l-3} t \nabla m. \end{aligned} $$

This accounts for (36).

For $i = k$,

$$ (47) \qquad \begin{aligned} \mathbb{E}(V(1:k)\bar{V}(r-(l-k)+1:r)V(1:k)) \\ &= \mathbb{E}(V(1:k)^2)\mathbb{E}(\bar{V}(r-(l-k)+1:r)) \\ &= \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{k-1} s^2 \left(\frac{m}{x_0} - 1\right)^{l-k} \\ &= \rho^{k-1} s^2 \left(\frac{m}{x_0} - 1\right)^{l+k-2}. \end{aligned} $$

This accounts for (37).
For $i > k$,

$$
\begin{align*}
& \mathbb{E}(V(1:k)\bar{V}(r-(l-i)+1:r)V(1:i)) \\
&= \mathbb{E}(V(1:k)\bar{V}(1:k))\mathbb{E}(V(k+1:i))\mathbb{E}(\bar{V}(r-(l-i)+1:r)) \\
&= \left(\frac{\sigma^2}{x_0^2} + \left(\frac{m}{x_0} - 1\right)^2\right)^{k-1} t \left(\frac{m}{x_0} - 1\right)^{i-k-1} \nabla m \left(\frac{m}{x_0} - 1\right)^{l-i} \\
(48) \qquad & = \rho^{k-1} \left(\frac{m}{x_0} - 1\right)^{l+k-3} t \nabla m.
\end{align*}
$$

This accounts for (38).

The constituent terms of all the separate sums of $\mathbb{E}(W_{r,k}^C W_{r,l}^C)$ stated in the Lemma are given by (44), (45), (43), (42), (46), (47) and (48). $\square$

The second moments can be bounded following a similar argument to that in Lemma 13. We obtain a general bound for the cycling estimate which does not require any assumption besides Assumption 1 and finiteness of the second moments; however, under Assumption 3, $\beta_0 < 1$ and the bound below could be simplified.

**Lemma 19.** For any $r \ge l \ge k$, $\text{cov}[W_{r,k}^C W_{r,l}^C]$ satisfies

$$
\begin{align*}
\operatorname{cov}[W_{r,k}^C W_{r,l}^C] &= \mathbb{E}(W_{r,k}^C W_{r,l}^C) - \mathbb{E}(W_{r,k}^C) \mathbb{E}(W_{r,l}^C) \\
&\le \frac{2}{r} \rho(\theta)^{k-1} \beta_0^{k+l-3} \max(\beta_0, 1) s(\theta)^2 (2l - 2).
\end{align*}
$$

*Proof.* We begin by bounding all the terms in the expression for $\mathbb{E}(W_{r,k}^C W_{r,l}^C)$ given in Lemma 18. Using the Cauchy-Schwarz inequality, $|t(\theta)| \le s(\theta)\beta_0 < s(\theta)$. By Jensen's inequality, $\nabla m(\theta)^2 \le s(\theta)^2$ and thus $|t(\theta)\nabla m(\theta)| < s(\theta)^2$.

$$
\text{term (33)} \le \frac{1}{r} \rho(\theta)^{k-1} \beta_0^{k+l-3} s(\theta)^2 \max(\beta_0, 1),
$$

$$
\text{term (34)} \le \frac{1}{r} \rho(\theta)^{k-2} \beta_0^{k+l-3} s(\theta)^2 (k-2),
$$

$$
\text{term (36)} \le \frac{1}{r} \rho(\theta)^{k-2} \beta_0^{k+l-3} s(\theta)^2 (k-1),
$$

$$
\text{term (37)} \le \frac{1}{r} \rho(\theta)^{k-1} \beta_0^{k+l-3} s(\theta)^2 \max(\beta_0, 1),
$$

$$
\text{term (38)} \le \frac{1}{r} \rho(\theta)^{k-1} \beta_0^{k+l-3} s(\theta)^2 (l - k - 1).
$$

Note that

$$
\begin{align*}
-\mathbb{E}(W_{r,k}^C) \mathbb{E}(W_{r,l}^C) &= -\left(\frac{m(\theta)}{x_0} - 1\right)^{k+l-2} (\nabla m(\theta))^2 \\
&\le \beta_0^{k+l-3} \max(\beta_0, 1) s(\theta)^2.
\end{align*}
$$

Now subtract $\mathbb{E}(W_{r,k}^C) \mathbb{E}(W_{r,l}^C)$ from $\mathbb{E}(W_{r,k}^C W_{r,l}^C)$.

$$
\text{term (35)} - \frac{1}{r} \sum_{i=1, i \ge k}^{r-l+1} \mathbb{E}(W_{r,k}^C) \mathbb{E}(W_{r,l}^C) = 0.
$$
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Comparison of work-normalised variance in a low, moderate and high variance scenario for several values of $n_0$. The stars denote the value of $\mathbb{E}[R]$ minimising the upper-bound of the WNV and the vertical line $\mathbb{E}[R] = n_0$.</figcaption></figure>

All the remaining (non-zero) terms of the subtraction of $\mathbb{E}(W_{r,k}^C) \mathbb{E}(W_{r,l}^C)$ from $\mathbb{E}(W_{r,k}^C W_{r,l}^C)$ may be bounded by summing the bounds above for (33), (34), (36), (37), (38) and $-\mathbb{E}(W_{r,k}^C) \mathbb{E}(W_{r,l}^C) (l+k-2)/r$ to obtain

$$ \mathbb{E}(W_{r,k}^C W_{r,l}^C) - \mathbb{E}(W_{r,k}^C) \mathbb{E}(W_{r,l}^C) \le \frac{2}{r} \rho(\theta)^{k-1} \beta_0^{k+l-3} \max(\beta_0, 1) s(\theta)^2 (l+k-2). $$

Recalling that $l \ge k$, we have the result. $\square$

<figure><img src="image_8.png" /></figure>

## APPENDIX C. ADAPTIVE TUNING OF p

As discussed in Section 3.2, the choice of $p$ is crucial in obtaining estimators which do not exhibit high variance. However, the results in Section 2 only give the range in which $p$ should be to guarantee finite variance, $(0, 1 - \beta^2)$.

Ideally, one would like to select $p$ that minimises the work-normalised variance $\text{WNV} = (n_0 + \mathbb{E}[R]) \times \text{var}[\hat{f}]$. However, in general, $\text{var}[\hat{f}]$ is unknown. One could replace $\text{var}[\hat{f}]$ with the upper-bounds established in Section 2 and minimise instead an upper-bound on WNV (Blanchet and Glynn, 2015). The upper bounds for the
simple and cycling estimators are

$$ 
\begin{aligned}
& (n_0 + \mathbb{E}[R]) \times \mathrm{var}[\hat{f}^S] \\
& \leq c^2(n_0 + (1-p)/p) \left[ \frac{\beta_0^2}{(1-\beta_0)^2} \times \frac{p}{1-p-\beta_0^2} + \frac{(1+\beta)}{1-\beta} \times \frac{(1-p)}{1-p-\beta^2} \right] \\
& (n_0 + \mathbb{E}[R]) \times \mathrm{var}[\hat{f}^C] \\
& \leq c^2(n_0 + (1-p)/p) \left[ \frac{\beta_0^2}{(1-\beta_0)^2} \times \frac{p}{1-p-\beta_0^2} + \frac{4p \log(1/p)}{(1-p-\beta^2)^2} \left[ \beta^2 + \frac{2\beta_0}{1-\beta_0^2} \right] \right]
\end{aligned}
$$

and are obtained by combining Proposition 4 and 7 and Proposition 4 and 9, respectively. In the interval $(0, 1 - \beta^2)$, both bounds are convex functions of $p$ whose unique minimum can be found using standard optimisation routines.

We compare the optimal value of $p$ obtained by minimising the bounds above in three scenarios: we set the mean of the $X_i$'s to be $m = 1$ and consider $\sigma^2 = 0.1, 1, 10$ which correspond to a low, moderate and high variance scenario. We consider three values of $n_0$ for the pilot run. Figure 6 shows the results for the simple and the cycling estimators. As expected, higher input variance requires higher $\mathbb{E}[R]$ to control the variance of the estimators. Generally, the optimal $p$ is larger for the cycling estimator, this is a consequence of the slower rate of divergence shown in Proposition 9.

We can see that, in the low and moderate variance case, our strategy of setting $\mathbb{E}[R] = n_0$ corresponds to selecting $p$ close to the optimal value for the cycling estimator.

## APPENDIX D. MVUE

Let $Y_i := X_i/x_0 - 1$, then (13) may be rewritten as $U_{r,k}^{\text{MVUE}} = S_{r,k}/(r_k)$, with

$$ S_{r,k} := \sum_{1 \le i_1 < \dots < i_k \le r} Y_{i_1} \dots Y_{i_k}, $$

for $r \ge k \ge 1$, which may be computed recursively as follows:

$$ S_{r,k} = S_{r-1,k} + Y_r S_{r-1,k-1} $$

with initialisations $S_{r,k} := 0$ for $r < k$, and $S_{r,0} := 1$ for $r \ge 0$. For a fixed $r$, it is thus possible to compute $U_{r,k}^{\text{MVUE}}$ for $k = 1, \dots, r$ at a $O(r^2)$ cost. This is the same computational complexity as the cycling method.

Figure 7 compares the variability of the simple, cycling and MVUE estimators in the toy model example of Section 4.1, when $d = 2$, $n = 1$ and $p = 10^{-2}$. One sees that cycling and MVUE seems to lead to the same level of performance.

## APPENDIX E. ADDITIONAL DETAILS FOR THE NUMERICAL EXPERIMENTS

E.1. **Toy LVM: MLMC estimators.** For a given $y \in \{y_1, \dots, y_n\}$ and a given $\theta$, the exact expressions for the moments of $X_j = p(y|Z_j, \theta)$ and $\beta^2$ are

$$ 
\begin{aligned}
m &= p(y|\theta) = \mathcal{N}(y; \theta \mathbf{1}_d, 2|\mathbf{d}_d), \\
\sigma^2 &= (4\pi)^{-d/2} \mathcal{N}(y; \theta \mathbf{1}_d, \frac{3}{2}|\mathbf{d}_d) - \mathcal{N}(y; \theta \mathbf{1}_d, 2|\mathbf{d}_d)^2 \\
\beta^2 &= 1 - \left(\frac{4}{3}\right)^{-d/2} \exp\left(-\frac{\|y - \theta \mathbf{1}_d\|^2}{6}\right).
\end{aligned}
$$
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Box-plots of 10³ independent realisations of the simple, cycling, and MVUE estimators for the toy LVM model when $d = 2$, $n = 1$ and $p = 10^{-2}$. The true value is indicated by the red dashed horizontal line.</figcaption></figure>

The maximum likelihood estimator is given by $\theta^* = (dn)^{-1} \sum_{i=1}^{n} 1_d y_i$.

Under the assumption that $\tilde{R} \sim \text{Geometric}(p)$, the MLMC estimators have finite variance and finite expected cost whenever $0.5 < p < 3/4$ (Shi and Cornish, 2021, Theorem 2 and 3). However the second moment of the cost is infinite. To see this consider the second moment of the cost

$$ \sum_{k=0}^{\infty} 2^{2(j+1+k)} \mathbb{P}(\tilde{R} = k) = \sum_{k=0}^{\infty} 2^{2(j+1+k)} p(1-p)^k = \frac{4^{j+1}p}{1-4(1-p)}, \quad \text{provided } p > 3/4. $$

Hence, it is not possible to guarantee that MLMC estimators have finite expectation and cost with finite variance at the same time.

E.2. **Toy LVM: High dimensions.** Figure 8 shows the results obtained with simple and cycling estimators, MLMC in the case $d = 20$ with $p = 1/97$ which gives $\mathbb{E}[R] = 96$. The cost of the simple and cycling estimators is stable but the estimate itself is not.

To further probe the cause of the large variance of the Taylor-based sum estimators, we separate the contribution to the variance due to the random truncation of the Taylor series, given by (7), from that due to the variance of the inputs $X_1, \dots, X_r$ given by (9). To find this contribution to the variance, we replace the inputs $X_j = p(y|Z_j, \theta)$ with inputs with zero variance, i.e $X_j \equiv p(y|\theta)$. In this case, simple and cycling estimators coincide and (9) is zero.

As can be seen from Figure 8, the estimates obtained when $X_j = p(y|\theta)$ overlap with the estimates with $X_j = p(y|Z_j, \theta)$, which implies that the poor performance of the Taylor-based sum estimators is attributable as much to the random truncation of the Taylor series as it is to the large variance of the $X$'s. In particular, as $d$ increases, the expansion point $x_0$ is no-longer close to $m$, and so $\beta^2$ approaches 1. With our choice of $p$, (7) is not finite. In order to guarantee finite variance for all $i = 1, \dots, n$ we would need to set $\mathbb{E}[R]$ of the order $10^8$.
<figure><img src="image_10.png" /><figcaption>FIGURE 8. Estimates of the marginal log-likelihood for a full data set of size $n = 1000$ plotted against computational cost (measured in number of samples $Z_i$ drawn) for an expected cost $E[R] = 96$ samples per data point and dimension $d = 20$. The vertical dashed line denotes the expected cost while the horizontal one denotes the true value of $\sum_{i=1}^{n} \log p(y_i|\theta)$.</figcaption></figure>

Similar considerations apply to MLMC, in fact, for this example we have KL($p(\cdot|y, \theta)$, $q(\cdot; y)$) = $\frac{d}{2}(\log \frac{4}{3} - \frac{1}{4}) \to \infty$ as $d \to \infty$, hence the importance weights get more and more degenerate as $d$ increases.

**E.3. Independent component analysis.** Recall that $\theta = (A, \sigma^2)$,

$$p(z|\theta) = p(z) = \prod_{j=1}^{d_z} \{\alpha\text{Logistic}(z_j; 1/2) + (1-\alpha)\delta_0(z_j)\},$$ $$p(y|z, \theta) = \mathcal{N}(y; Az, \sigma^2\mathbf{1}_{d_y}).$$

Given that the prior does not depend on $\theta$, we can write

$$\nabla_{\theta} p(y|\theta) = \int \nabla_{\theta} \log p(y|z, \theta) p(y|z, \theta) p(z) dz,$$

where $\nabla_{\theta} \log p(y|z, \theta)$ may be decomposed into

$$\nabla_{\sigma} \log p(y|z, \theta) = \frac{1}{\sigma^3} \|y - Az\|^2 - \frac{d_y}{\sigma}, \quad \nabla_{A} \log p(y|z, \theta) = \frac{1}{\sigma^2} (y - Az)z^T.$$

We select the following proposal distribution, as it is proportional to the likelihood $p(y|z, \theta)$:

$$(49) \qquad q(z; y, \theta) = \mathcal{N}(z; (A^T A)^{-1} A^T y, \sigma^2 (A^T A)^{-1}),$$
where again $\theta = (A, \sigma^2)$. The corresponding importance weights equal

$$w(z; y, \theta) := \frac{(2\pi\sigma^2)^{(d_z-d_y)/2}}{\det(A^T A)^{1/2}} \exp\left(-\frac{\|y\|^2 - y^T A(A^T A)^{-1} A^T y}{2\sigma^2}\right) p(z).$$

For the Taylor-based estimators, the $X_i$'s are obtained by sampling $Z_i$ from $q$ in (49) and setting $X_i = w(Z_i; y, \theta)$. To obtain the unbiased estimators of $\nabla_{\theta}p(y|\theta)$ we set $G_i = w(Z_i; y, \theta)\nabla_{\theta}\log p(y|Z_i, \theta)$, with $Z_i$ sampled from $q$.

An alternative and simpler proposal is given by $q(z) = p(z)$. However, this proposal leads to $(X_i, G_i)$'s with high variance which require large $p$ to ensure finite variance of the Taylor-based gradient estimators. On the contrary, with (49) it is sufficient to take $p \in (0, 0.9)$ to guarantee finite variance, and thus our choice of $p = 1/n_0 = 0.1$ is valid. Similar considerations apply to MLMC with proposal $q(z) = p(z)$, the estimates of the gradient are unstable due to the high variability of the weights.

For SAEM, we draw the initial value of $\sigma$ from a Uniform distribution over the interval $[0, 1]$ and the initial distribution of $A = (a_1, a_2)$ from a Uniform distribution over $[0, 1]^{512}$. We select the learning rate for SAEM to be $\Delta_t = 1/t$, with $t$ denoting the iteration number, which satisfies the conditions in Delyon et al. (1999) to guarantee convergence. The MCMC algorithm targeting the posterior $p(z|y, \theta)$ is the Metropolis-within-Gibbs algorithm described in Allassonnière and Younes (2012, Section 3.3).

Given the high dimension of the parameter space, gradient-based methods struggle if the starting point of the iteration is far from the optimum. Hence, we initialise all SGD algorithms at the value of $\theta$ obtained after one iteration of SAEM (Polyak and Juditsky, 1992). We use two different learning rates for $\sigma$ and $A$, since the gradients of the two parameters are of different order of magnitude. We set $\Delta_{\sigma} = 10^{-7}$ and $\Delta_{A} = 10^{-4}$. This is in line with the suggestion of Beatson and Adams (2019, Section 5.1). We stop all algorithms after 5000 iterations to guarantee convergence of all algorithms.

At each iteration of SGD, we set $x_0$ using $n_0 = 10$ pilot runs and obtain the corresponding $\beta^2$ as described in Section 3.2. If $1/n_0 \in (0, 1 - \beta^2)$, we set $p = 1/n_0$, i.e. $E[R] = n_0$; otherwise we set $p = (1 - \beta^2)/2$. For MLMC we adjust $j, p$ so that $E[2^{j+1+\tilde{R}}] \approx n_0$.

For SGD we use mini-batches of size 2, while for SAEM we scan the whole data set and, as suggested in Allassonnière and Younes (2012, Section 3.3), randomly update one component of the latent vector $z$ for each data point $y$. This can be seen as a form of mini-batch sampling as noticed in Kuhn et al. (2020).

The MSE for the parameters $(\sigma, A)$ is defined as $MSE(\sigma) := E[(\hat{\sigma} - \sigma)^2]$,

$$MSE(A) := E\left[\frac{1}{d_y}\sum_{i,j=1}^{d_y^{1/2}}((\hat{a}_1 - a_1)_{ij})^2 + \frac{1}{d_y}\sum_{i,j=1}^{d_y^{1/2}}((\hat{a}_2 - a_2)_{ij})^2\right],$$

where $\hat{\sigma}, \hat{a}_1, \hat{a}_2$ are the estimates obtained with each method.

To assess convergence of SAEM and SGD we check the value of $\sigma$ and of the trace of $a_1$ and $a_2$ over iterations (Figure 9).

**E.4. Exponential random graphs.** For a given $\theta$, we obtain an unbiased estimate of $\mathcal{Z}(\theta)$ by running a waste-free SMC sampler (Dau and Chopin, 2022) based on adaptive tempering, i.e. the intermediate distributions are $\pi_t(y) \propto \exp\{\lambda_t s(y)\}$,
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Converge of SAEM and SGD for the ICA model.</figcaption></figure>

with $0 = \lambda_0 < \dots < \lambda_T = \theta$, and at beginning of iteration $t$, $\lambda_t$ is set so that ESS (effective sample size) equals $N/2$, where $N$ is the number of particles. Since we use a waste-free sampler, $N = M \times P$, where $M$ is the number of resampled particles, here $M = 100$, and $P$ is the length of the simulated chains, here $P = 500$. For the low variance setting, we use the same approach except $M = 10^3$. Note that the initial distribution $\pi_0$ is easy to sample from, as it is the uniform distribution over the $\binom{n}{2}$ possible networks. The Markov kernel used at time $t$ to move the particles is a Metropolis kernel based on the following proposal: one edge $(i, j)$ is chosen at random, and $y_{ij}$ is flipped. We use the Python package particles (https://github.com/nchopin/particles) to implement this SMC sampler.

The proposal $q$ is Gaussian with mean $\mu = (-3.27, 0.6)$ and covariance matrix

$$ \Sigma = \begin{pmatrix} 0.2016 & -0.0432 \\ -0.0432 & 0.0108 \end{pmatrix}. $$

This proposal was obtained through a small number of pilot runs.

The $\theta_i$'s are generated using a scrambled Sobol' sequence, see Owen (1997).

### E.5. Exponential random graphs: Comparison with the exchange algorithm.
We compare our approach with one of the state of the art for Bayesian inference on ERG models, i.e. the exchange algorithm of Caimo and Friel (2011) implemented in the bergm package in R. We use the same setup described in Section 4.5 therein. The shape of the posterior distribution seems to be in good agreement with the obtained in our experiments (see Figure 10).

As we now discuss, there are good reasons in favour of our method. First, the exchange algorithm of Caimo and Friel (2011) requires running nested MCMC chains where the inner one samples from the likelihood $p(y|\theta)$; in addition, in the version of the algorithm implemented in bergm, population Monte Carlo is employed to improve mixing. As a consequence, the resulting method is rather expensive and, due to the serial nature of MCMC, harder to parallelise than our proposed approach, in which each $\mathcal{Z}(\theta)$ is obtained from different independent runs of SMC which can run in parallel.

Secondly, a fundamental task in the context of network model is the ability to select the model that best fits the data. This type of inference requires estimation of the model evidence $p(y)$. To the best of our knowledge, given the output of the
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Bivariate weighted histograms approximating the posterior distributions obtained with the simple estimator, the cycling estimator and the bergm package implementing the exchange algorithm of Caimo and Friel (2011).</figcaption></figure>

exchange algorithm one can compute the model evidence

$$p(y) = p(\theta) \frac{\exp\{\theta^T s(y)\}}{\mathcal{Z}(\theta)} \frac{1}{p(\theta|y)}$$

in two ways:

(1) An approach similar to ours in which $\log \mathcal{Z}(\theta)$ is estimated via path sampling, and $p(\theta|y)$ via the exchange algorithm and a kernel density estimate (KDE) (Caimo and Friel, 2013, Section 4).

(2) Pseudo-likelihood approaches in which the intractable likelihood is replaced by a tractable approximation (Bouranis et al., 2018).

The first approach requires running the exchange algorithm to sample from $p(\theta|y)$, from which a KDE is formed. Path sampling with auxiliary MCMC chains is used to sample from $p(y|\theta)$. This approach is clearly only suitable for a small number of parameters due to the instability of the KDE in high dimensions. Also, as mentioned by the authors, this strategy incurs two sources of error when estimating $\log \mathcal{Z}(\theta)$: the Monte Carlo error to approximate the expectation w.r.t. $p(y|\theta)$, and the error due to discretising the path integral. The resulting approximation of the model evidence is naturally biased. Our approach allows computation of the model evidence directly from the samples used to estimate the posterior, without requiring additional samples from $p(y|\theta)$ and without the use of KDE. It also provides unbiased estimators which can be used within pseudo-marginal methods for model selection.

The second approach also cannot provide unbiased estimators of the model evidence and the quality of the approximation largely depends on the pseudo-likelihood used. To obtain a likelihood approximation close to the true one, Bouranis et al. (2018) propose several improvements. We tested the evidence estimation in the bergm package on our toy example and found that the Monte Carlo MLE estimator needed to implement their proposed adjustments to the pseudo-likelihood could not converge since the underlying MCMC was not mixing.

E.6. **Comparison with biased estimators.** To further motivate the need for unbiased estimators, we compare the MLMC estimator of Shi and Cornish (2021) and our simple and cycling estimators with their biased counterparts.
<figure><img src="image_13.png" /><figcaption>FIGURE 11. Mean absolute deviation and bias over 100 repetitions against cost for 2 biased estimators (IWAE, naive) and 3 unbiased estimators (MLMC, simple, cycling). The vertical dashed lines are the expected costs used in Figure 2 of the revised/original manuscript, namely $E[R] = 6$ and $E[R] = 96$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12. Distribution of the estimates of $f(m)$ for small ($E[R] = 6$) and high ($E[R] = 96$) cost. The dashed horizontal line is the true value.</figcaption></figure>

E.6.1. *Toy LVM.* We first consider the toy LVM and two biased estimators:

* the naive estimator $f(\bar{X})$ where $\bar{X}$ is the mean of the $X_i$'s;

* the IWAE estimator $I_j$ given below (18);

and their unbiased counterparts: simple/cycling and MLMC, respectively (see Figure 11 and Figure 12).

Figure 11 shows the behaviour of the relative mean absolute deviation $\mathbb{E}[|f(m) - \hat{f}|]/|f(m)|$ and bias as the expected computational cost increases for the toy LVM model with $d = 2$.

For large enough cost all estimators except the simple one have very small relative mean absolute deviation. This is because the variance of the simple estimator does not converge to zero as the cost increases, as shown in Example 8. However, for small cost the average (out of 100 replicates) bias of both IWAE and the naive estimator is considerably larger than that of the cycling estimator (see Figure 11).

Figure 12 shows the distribution of the 100 replicates for the smallest cost (left-most vertical dashed line) and the highest one. We can see that cycling provides a provably unbiased estimator without increasing the variability in the original biased naive estimator. When the cost is large it is also clear that cycling outperforms simple since the variability of the latter does not decrease to zero as cost increases.
<table><caption>TABLE 3. Comparison of reconstruction accuracy and cost for the censored logistic ICA model. The results are averaged over 100 repetitions.</caption><thead><tr><th>Method</th><th>MSE(σ)</th><th>MSE(A)</th><th>runtime (s)</th></tr></thead><tbody><tr><td>SAEM</td><td>7.00 · 10<sup>-7</sup></td><td>0.22</td><td>33</td></tr><tr><td>MLMC</td><td>2.95 · 10<sup>-6</sup></td><td>0.29</td><td>11</td></tr><tr><td>simple</td><td>2.97 · 10<sup>-6</sup></td><td>0.27</td><td>10</td></tr><tr><td>cycling</td><td>6.80 · 10<sup>-5</sup></td><td>0.24</td><td>14</td></tr><tr><td>IWAE</td><td>3.55 · 10<sup>-6</sup></td><td>0.30</td><td>7</td></tr><tr><td>Plug in</td><td>3.62 · 10<sup>-6</sup></td><td>0.27</td><td>7</td></tr></tbody></table>

<table><caption>TABLE 4. Work-normalised variance and mean squared error for $10^3$ replicates of $\tilde{Z}(\theta)^{-1}$ for the simple, the cycling estimator and the plug in estimator.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Simple</th><th colspan="2">Cycling</th><th colspan="2">Plug in</th></tr><tr><th>WNV</th><th>MSE</th><th>WNV</th><th>MSE</th><th>WNV</th><th>MSE</th></tr></thead><tbody><tr><td>Moderate</td><td>61</td><td>0.32</td><td>24</td><td>0.15</td><td>26</td><td>0.29</td></tr><tr><td>Low</td><td>83</td><td>0.04</td><td>26</td><td>0.01</td><td>105</td><td>0.09</td></tr></tbody></table>

### E.6.2. Independent component analysis.
We also compare the plug in estimator $f(\bar{X})$ and the IWAE one for the ICA model using the same setup of Table 3 where we set the number of samples for IWAE and the plug in estimator to 10.

### E.6.3. Exponential random graphs.
To further showcase the strengths of the cycling estimator we also consider the ERG model. In particular, the setup of Table 2 for $E[R] = 10$. We compare the work-normalised-variance and the relative mean squared error (MSE) for simple, cycling and the plug in estimator with same number of replicates $R$. The relative MSE is obtained by considering as true value of $1/Z(\theta)$ that obtained by taking the reciprocal of 30,000 unbiased values of $Z(\theta)$ obtained in the low variance regime (see Table 4). It is evident that in both a low variance regime and in a medium one cycling is preferable to the simple plug in estimator. In the low variance regime even the simple estimator is preferable to the plug in one.

ENSAE, INSTITUT POLYTECHNIQUE DE PARIS, FRANCE

Email address: nicolas.chopin@ensae.fr

ESOMAS, UNIVERSITÀ DI TORINO, ITALY

Email address: francescaromana.crucinio@unito.it

NIASRA AND SCHOOL OF MATHEMATICS AND APPLIED STATISTICS, UNIVERSITY OF WOLLONGGONG, AUSTRALIA

Email address: sumeetpals@uow.edu.au