A key challenge in applying Monte Carlo and variational inference (VI) is the design of proposals and variational families that are flexible enough to closely approximate the posterior, but simple enough to admit tractable densities and variational bounds.
This paper presents recursive auxiliary-variable inference (RAVI), a new framework for exploiting flexible proposals, for example based on involved simulations or stochastic optimization, within Monte Carlo and VI algorithms.
The key idea is to estimate intractable proposal densities via meta-inference: additional Monte Carlo or variational inference targeting the proposal, rather than the model.
RAVI generalizes and unifies several existing methods for
inference with expressive approximating families,
which we show correspond to specific choices
of meta-inference algorithm, and provides new theory for analyzing their bias and variance.

We illustrate RAVI's design framework and theorems by using them to analyze and improve upon Salimans et al., 2015's Markov Chain Variational Inference, and to design a novel sampler for Dirichlet process mixtures, achieving state-of-the-art results on a standard benchmark dataset from astronomy and on a challenging data-cleaning task with Medicare hospital data.

§ INTRODUCTION

Monte Carlo and variational inference algorithms
are the workhorses of modern probabilistic inference,
a fundamental problem with applications in many disciplines [Murphy, 2012].
A key challenge in applying these algorithms is the design of proposal distributions (in VI, variational families), which can greatly affect their performance [Chatterjee et al., 2018]. A good proposal should incorporate any knowledge the practitioner might have about the shape of the posterior; however, this goal is often in tension with the requirement that a proposal's marginal density be analytically tractable, in order to compute importance weights, MCMC acceptance probabilities, or gradient updates for VI. The challenge is that proposal distributions that are simple enough to admit tractable densities may not be flexible enough to solve real-world posterior inference problems.

In this paper, we present a new framework, called Recursive Auxiliary-Variable Inference (RAVI), for incorporating more complex proposals, without tractable densities, into standard Monte Carlo and VI algorithms.
The key idea is to approximate the intractable
proposal densities using meta-inference [Cusumano-Towner and Mansinghka, 2017]: an additional layer of Monte Carlo or variational inference targeting the proposal, rather than the model. RAVI generalizes and unifies several existing methods for inference with expressive proposals [Salimans et al., 2015, Ranganath et al., 2016, Sobolev and Vetrov, 2019], which we show correspond to specific choices of meta-inference algorithm (see Table <ref> and Appendix <ref> for 13 examples).

Contributions. Our key contributions are:

* the RAVI framework, including new recursive algorithms for IS, VI, SMC, and MH using proposals that lack tractable marginal densities (Sections <ref> and <ref>);

* theorems that precisely characterize the impact of RAVI's approximations on overall inference quality (sampler variance, or tightness of the variational bound) (Section <ref>); and

* two extended examples of RAVI's application to algorithm design and analysis: (1) a novel variant of Salimans et al., 2015's Markov Chain Variational Inference (MCVI) algorithm that, unlike vanilla MCVI, scales to handle proposals incorporating long MCMC chains; and (2) a novel sampler for Dirichlet process mixtures that uses a randomized agglomerative clustering algorithm as a proposal, outperforming strong baselines on a standard benchmark from astronomy [Drinkwater et al., 2004] and a challenging data-cleaning problem using Medicare data [Medicare, 2012], adapted from Lew et al., 2021.

§ RECURSIVE AUXILIARY-VARIABLE INFERENCE

Structure of a RAVI inference strategy $\mathcal{S}$ targeting the posterior $p(x \mid y)$ of a latent-variable model.
The proposal $\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x) \text{d}r$ has an intractable marginal density, so the strategy also specifies a meta-inference strategy $\mathcal{S}.\mathcal{M}$ that targets $\mathcal{S}.q(r \mid x)$. Nesting continues until the $q$
approximation at some layer has a tractable density, at which
point no further meta-inference is needed.

Intractable posterior approximation $\mathcal{S}.q$
Meta-posterior approximation $\mathcal{S}.\mathcal{M}(x).q$

IWAE [Burda et al., 2015]
t]@l@SIR with tractable proposal
$q(x) = \sum_{j=1}^K \int \prod_{i=1}^K q(x_i) \frac{w_j}{\sum_{l=1}^K w_l}\delta_{x_j}(x)\text{d}x_{1:K}$     t]@l@Conditional SIR
$ h(x_{1:K}, j) = \frac{1}{K} \delta_{x}(x_j) \prod_{i \neq j} q(x_i)$

HVM [Ranganath et al., 2016]       t]@l@Latent-variable proposal
$q(x) = \int q(z, x)\text{d}z$
t]@l@Tractable proposal
$h(z; x)$

MCVI [Salimans et al., 2015]     t]@l@MCMC chain initialized from tractable proposal
$q(x) = \int q(x_0) \prod_{i=1}^M T(x_i; x_{i-1}) \delta_{x_M}(x) \text{d}x_{0:M}$
t]@l@Sequence of tractable proposals
$h(x_{0:M}; x) = \delta_{x}(x_M) \prod_{i=0}^{M-1} r_i(x_i; x_{i+1})$

IWHVI [Sobolev and Vetrov, 2019]     t]@l@Latent-variable proposal
$q(x) = \int q(z, x) \text{d}z$
t]@l@SIR with tractable proposal

$h(z; x) = \sum_{j=1}^K \int \prod_{i=1}^K h(z_i; x) \frac{w_j}{\sum_{l=1}^K w_l}\delta_{z_j}(z)\text{d}z_{1:K}$

DIWHVI [Sobolev and Vetrov, 2019]     t]@l@SIR with latent-variable proposal
$q(x) = \sum_{j=1}^K \int \prod_{i=1}^K q(z_i, x_i)\cdot$
$\frac{w_j}{\sum_{l=1}^K w_l} \delta_{x_j}(x) \text{d}(z, x)_{1:K}$
t]@l@Conditional SIR for $j, x_{1:K}, z_{-j}$

$h(z_{-j}, x_{1:K}, j; x) = \frac{1}{K} \delta_{x}(x_j) \prod_{i \neq j} q(z_i, x_i)$

SIR with tractable proposal for $z_j$

$h(z_j; x) = \sum_{k=1}^K \int \prod_{i=1}^K h(z_j^{(i)}; x) \frac{w_k}{\sum_{m=1}^K w_m}\delta_{z_j^{(k)}}(z_j)\text{d}z_j^{(1:K)}$

Many existing approaches to inference with flexible proposals can be understood as
instances of RAVI, with particular choices of inference strategy $\mathcal{S}$.
In some cases, including IWHVI and DIWHVI, the meta-posterior approximation is itself intractable,
necessitating meta-meta-inference $\mathcal{S}.\mathcal{M}(x).\mathcal{M}$. See Appendix <ref> for eight additional examples.

In this section, we introduce the RAVI framework in the context of a running example:
we incorporate a chain of MCMC steps into a proposal, so that it can more accurately approximate a
posterior distribution. Our approach generalizes Salimans et al., 2015's Markov Chain Variational Inference (MCVI) algorithm, and fixes a flaw that prevents it from scaling to longer MCMC chains.

An expressive proposal with an intractable marginal density.
Let $p(x, y)$ be a latent-variable model and
$y$ an observation. We are interested in approximating
$p(x \mid y)$ using an expressive proposal $q(x)$, that
generates an initial location $x_0$ from a simple parametric distribution $q_0$, then
iterates $M$ steps of an MCMC kernel $T$:

$$q(x) = \int q_0(x_0) \prod_{i=1}^M T(x_{i-1}\rightarrow x_i) \delta_{x_M}(x) \text{d}x_{0:M}.$$

Even when $q_0$ is a poor approximation to $p(x \mid y)$, $q(x)$ may be close to the posterior, if $M$ is sufficiently high. However, because the density $q(x)$ cannot be tractably evaluated, we cannot use $q(x)$ as a proposal within importance sampling (we have no way to evaluate the importance weight $\frac{p(x, y)}{q(x)}$), nor as a variational family in VI (we cannot estimate the ELBO $\mathcal{L} = \mathbb{E}_{x \sim q} [\log \frac{p(x,y)}{q(x)}]$ or its gradient, making it impossible to learn $p$'s or $q$'s parameters).

Approximating intractable densities with meta-inference. RAVI's goal is to enable inference even when we cannot compute the marginal densities of our proposals and variational families. To apply RAVI, we must specify
not just the proposal itself
but a complete inference strategy:

Definition. An inference strategy $\mathcal{S}$ targeting $\pi$ specifies:

* a posterior approximation $\mathcal{S}.q(x) \approx \pi(x)$[To simplify the exposition, we assume that if an inference strategy $\mathcal{S}$ targets $\pi$, then the approximation $\mathcal{S}.q$ is mutually absolutely continuous with $\pi$, i.e. the measure-zero events under $\pi$ are exactly the same as those under $\mathcal{S}.q$. This requirement can be relaxed somewhat; see Appendix <ref>.] that either has a tractable density, or is the marginal distribution of a joint distribution with a tractable density, i.e. $\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x) \text{d}r$, and,

* if $\mathcal{S}.q$ has an intractable marginal density, a meta-inference strategy $\mathcal{S}.\mathcal{M}$, assigning to each value of $x$ an inference strategy $\mathcal{S}.\mathcal{M}(x)$ targeting $\mathcal{S}.q(r \mid x)$.

Figure <ref> illustrates the recursive
structure of an inference strategy.
The key novelty is the inclusion of meta-posterior approximations:
additional proposals that the user must specify for
inferring auxiliary variables introduced by existing proposal distributions. In our running example, we take $\mathcal{S}.q(x)$ to be our MCMC-based posterior approximation: it lacks a tractable density, but is the marginal of a tractable joint density $\mathcal{S}.q(x_{0:M}, x)$ over entire MCMC traces. A meta-posterior approximation, then, is a probability distribution $\mathcal{S}.\mathcal{M}(x).q(x_{0:M})$ that approximates the meta-posterior $\mathcal{S}.q(x_{0:M} \mid x)$: the distribution over traces of the MCMC chain, given
the final location $x$.

The meta-posterior approximations enable RAVI to estimate the intractable marginal density of the top-level posterior approximation, to compute weights and gradients:

In Monte Carlo: If $\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x) \text{d}r$ is intended for use as a Monte Carlo proposal, RAVI uses meta-inference to obtain an unbiased estimate of $\frac{1}{\mathcal{S}.q(x)}$ (Algorithm [alg:alg2]2), which is then multiplied by $p(x, y)$ to estimate the importance weight $\frac{p(x, y)}{\mathcal{S}.q(x)}$. This process relies on the harmonic mean identity [Newton and Raftery, 1994], that for any meta-posterior approximation $h$,
$$\mathbb{E}_{\mathcal{S}.q(r \mid x)}\left[\frac{h(r)}{\mathcal{S}.q(r, x)}\right] = \frac{1}{\mathcal{S}.q(x)}\mathbb{E}\left[\frac{h(r)}{\mathcal{S}.q(r \mid x)}\right] = \frac{1}{\mathcal{S}.q(x)}.$$

(Harmonic mean estimators are infamous for having potentially infinite variance, but only when $h$ is set to a broad prior; we give a general analysis of the variance of RAVI's importance weights in Section <ref>.)

In Variational Inference: If $\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x) \text{d}r$ is intended as a variational family, then RAVI uses the meta-posterior approximation to formulate an upper bound on $\log \mathcal{S}.q(x)$: for any meta-posterior approximation $h(r)$,

$$\log \mathcal{S}.q(x) \leq \mathcal{U}(x) := \mathbb{E}_{\mathcal{S}.q(r \mid x)}[\log \mathcal{S}.q(r, x) - \log h(r)].$$

This follows from Jensen's inequality, and the harmonic
mean identity from above.
With this upper bound in hand, we formulate a surrogate ELBO $\mathcal{L}_\mathcal{S} = \mathbb{E}_{\mathcal{S}.q(x)}[\log p(x, y) - \mathcal{U}(x)] \leq \mathcal{L}$, which we can tractably estimate and optimize via stochastic gradient descent (Algorithm [alg:alg3]3).

In Section <ref>, we show how similar estimators can be built up recursively when the meta-posterior approximations themselves have intractable marginal densities.

A meta-inference strategy that recovers the MCVI objective [Salimans et al., 2015]. In our running example, where the auxiliary randomness $r$ is a trace $x_{0:M}$ of locations visited by MCMC, one option for meta-inference is to learn neurally parameterized
reverse Markov kernels $R_i(x_{i+1} \rightarrow x_i)$, and apply them
in sequence to infer a plausible trace of MCMC steps leading to the
final location $x$:

$$\mathcal{S}.\mathcal{M}(x).q(x_{0:M}) = \delta_{x}(x_M) \prod_{i=0}^{M-1} R_i(x_{i+1} \rightarrow x_i).$$

This approximation to $\mathcal{S}.q(x_{0:M} \mid x)$ has a tractable
density, and so completely specifies the meta-inference strategy $\mathcal{S}.\mathcal{M}$; there is no need to specify a meta-meta-inference
strategy. Given $\mathcal{S}$, RAVI optimizes the surrogate objective
$\mathcal{L}_\mathcal{S} = \mathbb{E}_{x \sim \mathcal{S}.q}[\log p(x, y) - \mathcal{U}_{\mathcal{S}.\mathcal{M}(x)}]$, where

$$\mathcal{U}_{\mathcal{S}.\mathcal{M}(x)} = \mathbb{E}_{x_{0:M} \sim \mathcal{S}.q(x_{0:M} \mid x)}\left[\log \frac{\mathcal{S}.q(x_{0:M}, x)}{\mathcal{S}.\mathcal{M}(x).q(x_{0:M})}\right].$$

For the above choice of $\mathcal{S}.\mathcal{M}$, the RAVI objective $\mathcal{L}_\mathcal{S}$ exactly coincides with the Markov Chain Variational Inference (MCVI) objective of Salimans et al., 2015. In fact, RAVI unifies and generalizes many existing methods; 13 examples are collected in Table <ref> and Appendix <ref>.

Analyzing MCVI within the RAVI framework. Framing MCVI as a RAVI algorithm lets us analyze it using general theory about RAVI objectives. For example, the relative tightness of the bound $\mathcal{L}_\mathcal{S}$ is controlled by the quality of meta-inference:

$$\mathcal{L} - \mathcal{L}_\mathcal{S} = \mathbb{E}_{\mathcal{S}.q(x)}[KL(\mathcal{S}.q(x_{0:M} \mid x) || \mathcal{S}.\mathcal{M}(x).q(x_{0:M}))].$$

We can use this characterization to analyze the MCVI objective's behavior as $M$ grows, i.e., as MCMC steps are added. Informally, as the MCMC chain begins to mix, the marginal distribution $\mathcal{S}.q(x)$ over the final location of the chain should grow closer to the posterior $p(x \mid y)$, tightening the (intractable) ELBO $\mathcal{L}$. Unfortunately, the meta-inference gap $\mathcal{L} - \mathcal{L}_\mathcal{S}$ grows with $M$, unless each kernel $R_i$ exactly captures the local
posterior $\mathcal{S}.q(x_i \mid x_{i+1})$. (This can be seen as an instance of the well-known degeneracy problem of sequential importance sampling <cit.>.) As MCMC converges, the rate of improvement in $\mathcal{L}$ slows, and the meta-inference penalty
for increasing the chain's length eventually outweighs the
benefit of improving the posterior approximation $\mathcal{S}.q$. The red curves in Figure <ref> show this phenomenon playing out on two toy targets: we see that $\mathcal{L}_\mathcal{S}$ does
become tighter as more MCMC steps are added, but only to a point, before the bound begins to loosen.

Resolving the issue with improved meta-inference.
RAVI clarifies that the variational bound loosens with increasing $M$ due to poor meta-inference: as the MCMC chain grows longer, error in the learned backward kernels accumulates.
This analysis also points to
a solution: use a meta-inference
strategy $\mathcal{S}.\mathcal{M}$ that can scale to
longer MCMC histories.

RAVI Inference Strategy 1: RAVI-MCVI
mcvirmcvi($M, K$).qmcvimrmcvi($M, K$).M($x$).qmcvimmrmcvi($M, K$).M($x$).M($x_{0:M})$.q
infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

latent variable $x$
MCMC trace $x_{0:M}$
$x_0 \sim q_0$
$i \in 1, \dots, M$
$x_i \sim T(x_{i-1} \rightarrow \cdot)$

metaalgMeta-Posterior Approx.

MCMC trace $x_{0:M}$
SMC particles $x_{0:M}^{1:K}$, ancestor indices $a_0, a_{1:M}^{1:K}$
$k \in 1,\dots, K$
$(x_M^k, w_M^k, t_k) \gets (x, q_m(x), [x])$

$i \in M-1, \dots, 0$
$k \in 1,\dots,K$
$a_{i+1}^k \sim \text{Discrete}(w_{i+1}^{1:K})$
$x_i^k \sim R_i(x_{i+1}^{a_{i+1}^k} \rightarrow \cdot)$*MCVI backward kernel
$w_i^k \gets \frac{q_i(x_i^k)T(x_i \rightarrow x_{i+1}^{a_{i+1}^k})}{q_{i+1}(x_{i+1}^{a_{i+1}^k})R_i(x_{i+1}^{a_{i+1}^k} \rightarrow x_i^k)}$
$t_k \gets [x_i^k, t_k^{a_{i+1}^k}\dots]$

$a_0 \sim \text{Discrete}(w_{0}^{1:K})$

metaalgMeta-Meta-Posterior Approx.

SMC particles $x_{0:M}^{1:K}$, ancestor indices $a_0, a_{1:M}^{1:K}$

$i \in 0, \dots, M$
$b_i \sim  \text{Uniform}(1, K)$

$k \in 1,\dots, K$
$(x_M^k, w_M^k) \gets (x, q_m(x)])$

$i \in M-1, \dots, 0$
$k \in 1,\dots,K$
$k = b_i$
$(a_{i+1}^k, x_i^k) \gets (b_{i+1}, x_i)$

$a_{i+1}^k \sim \text{Discrete}(w_{i+1}^{1:K})$
$x_i^k \sim R_i(x_{i+1}^{a_{i+1}^k} \rightarrow \cdot)$

$w_i^k \gets \frac{q_i(x_i^k)T(x_i \rightarrow x_{i+1}^{a_{i+1}^k})}{q_{i+1}(x_{i+1}^{a_{i+1}^k})R_i(x_{i+1}^{a_{i+1}^k} \rightarrow x_i^k)}$

$a_0 \gets b_0$
$(a_0, a_{1:M}^{1:K}, x_{0:M}^{1:K})$

A standard approach
to resolving the degeneracy problem when inferring
sequences of latent variables is sequential
Monte Carlo (SMC) [Del Moral et al., 2006]. SMC tracks
$K$ hypotheses about a latent sequence, periodically
weighting the hypotheses and resampling, to
clone promising particles and cull poor ones. Using
RAVI, we can use SMC for
meta-inference: we choose $\mathcal{S}.\mathcal{M}(x).q(x_{0:M})$ to generate a collection of $K$ possible backward
MCMC trajectories, using SMC, before selecting one to return.
This meta-posterior approximation is shown in RAVI Inference Strategy [infstrat:mcvi]1.

This algorithm does not itself have a tractable marginal density: computing $\mathcal{S}.\mathcal{M}(x).q(x_{0:M})$

would require large sums over the resampling variables and intractable integrals over the particle collection. But
this is where RAVI's recursive structure comes into play:
a meta-inference strategy may use an intractable meta-posterior
approximation, so long as
we attach a meta-meta-inference strategy $\mathcal{S}.\mathcal{M}(x).\mathcal{M}(x_{0:M}).q(a_0, a_{1:M}^{1:K}, x_{0:M}^{1:K})$. In this case meta-meta-inference must infer the auxiliary variables of SMC (ancestor variables and unchosen trajectories), given the final chosen trajectory $x_{0:M}$. For this we can use the conditional SMC algorithm [Andrieu et al., 2010],
which runs SMC, with the same auxiliary variables, but constrained to ensure that one of the $K$ particles
traces the observed trajectory $x_{0:M}$. Because cSMC
introduces no new auxiliary variables, it has a
tractable density, and there is no need to specify
a fourth layer of meta-inference. The full
tower of posterior approximations is given in RAVI Inference Strategy [infstrat:mcvi]1.

In Section <ref>, we compare MCVI to , for different $K$ and $M$. Figure <ref> shows that meta-inference error is greatly reduced by using SMC, so that the variational bound $\mathcal{L}_\mathcal{S}$ continues to tighten as the MCMC chain grows longer.

Recursive Monte Carlo Estimation

Algorithm 1: RAVI Importance Sampling ()
unnormalized target $\tilde{\pi}(x) = Z\pi(x)$
inference strategy $\mathcal{S}$
$(x, \hat{Z})$ properly weighted for $\pi(x)$, s.t. $\mathbb{E}[\hat{Z}] = Z$
$\mathcal{S}.q$ has a tractable marginal density
$x \sim \mathcal{S}.q$
$w \gets \frac{1}{\mathcal{S}.q(x)}$

$\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x)\text{d}r$
$(r, x) \sim \mathcal{S}.q$
$w \gets \texttt{HME}(\mathcal{S}.q(\cdot \mid x), r, \mathcal{S}.\mathcal{M}(x))$

$(x, w\tilde{\pi}(x))$

Algorithm 2: RAVI Harmonic Mean Estimation ()
unnormalized target $\tilde{\pi}(x) = Z\pi(x)$
exact sample $x \sim \pi$
inference strategy $\mathcal{S}$
unbiased estimate $\check{Z}^{-1}$ of $Z^{-1}$
$\mathcal{S}.q$ has a tractable marginal density
$w \gets \mathcal{S}.q(x)$

$\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x)\text{d}r$
$(r, w) \gets \texttt{IMPORTANCE}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))$

Recursive Variational Objectives and Gradient Estimation

Algorithm 3: RAVI ELBO and gradient estimator ($\texttt{ELBO}\nabla$)
model $p(x, y)$
data $y$
inference strategy $\mathcal{S}$
unbiased estimates of $\mathcal{L}(p, y, \mathcal{S})$ and of $\nabla_\theta \mathcal{L}(p, y, \mathcal{S})$
$\mathcal{S}.q$ has a tractable marginal density
$x \sim \mathcal{S}.q$
$ (\hat{U}, \widehat{\grad}) \gets (\log \mathcal{S}.q(x), \grad \log \mathcal{S}.q(x) \cdot (1 + \log \mathcal{S}.q(x)))$
$\mathbf{g} \gets \grad \log \mathcal{S}.q(x)$

$\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x)\text{d}r$
$(r, x) \sim \mathcal{S}.q$
$(\hat{U}, \widehat{\nabla_\theta}, \mathbf{g}) \gets \texttt{EUBO}\nabla(\mathcal{S}.q, x, r, \mathcal{S}.\mathcal{M}(x))$

$\hat{L} \gets \log p(x, y) - \hat{U}$
$\widehat{\grad}'\gets \grad \log p(x, y)
+ \mathbf{g}
\log p(x, y)
- \widehat{\grad}.$
$(\hat{L}, \widehat{\grad}')$

Algorithm 4: RAVI EUBO and gradient estimator ($\texttt{EUBO}\nabla$)
model $p(x, y)$
data $y$
exact sample $x \sim p(x \mid y)$
inference strategy $\mathcal{S}$
unbiased estimates of $\mathcal{U}(p, y, \mathcal{S})$ and $\nabla_\theta \mathcal{U}(p, y, \mathcal{S})$
quantity $\mathbf{g}$ (see Thm. 2)
$\mathcal{S}.q$ has a tractable marginal density
$ (\hat{L}, \widehat{\grad}) \gets (\log \mathcal{S}.q(x), \grad \log \mathcal{S}.q(x))$

$\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x)\text{d}r$
$(\hat{L}, \widehat{\nabla_\theta}) \gets \texttt{ELBO}\nabla(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))$

$\hat{U} \gets \log p(x, y) - \hat{L}$
$\mathbf{g} \gets \nabla_\theta \log p(x, y)$
$\widehat{\grad}'\gets \grad \log p(x, y)
+ \mathbf{g} \cdot \hat{U}
- \widehat{\grad}$
$(\hat{U}, \widehat{\grad}', \mathbf{g})$

Using the inference strategy within a Monte Carlo algorithm, to estimate marginal likelihoods from MCMC results. Our inference strategy $\mathcal{S}$ can also be used as proposal within Monte Carlo algorithms, such as importance sampling. In the context of our example, where $\mathcal{S}.q$ incorporates $M$ steps of a Markov chain,
this allows us to assign an importance weight to
each run of the Markov chain.
The weight is an unbiased estimate of the marginal likelihood $p(y)$ of the model;

thus, we can view the algorithm as a way to derive marginal likelihood estimates from MCMC runs, a task of long-standing interest in the Monte Carlo community [Neal, 2001]. In Section <ref>, we show that in some settings MCVI compares favorably a standard algorithm for the task, annealed importance sampling (AIS) [Neal, 2001].

§ ALGORITHMS

In this section, we present algorithms for using RAVI inference strategies within Monte Carlo and variational inference algorithms, as proposals and variational families.

RAVI for Importance Sampling and SMC. In importance sampling and SMC algorithms, proposals $q$ are used to (1) generate proposed values $x \sim q$, and (2) compute importance weights $\frac{p(x)}{q(x)}$. But in both IS and SMC, it suffices to produce unbiased estimates of $\frac{p(x)}{q(x)}$ [Liu, 2008, Chopin and Papaspiliopoulos, 2020].

RAVI exploits this degree of freedom to generate proper importance weights even when $q(x)$ is intractable.
Suppose $\tilde{\pi} = Z\pi$ is an unnormalized target density, and $\mathcal{S}$ is a RAVI inference strategy targeting $\pi$. Algorithm [alg:alg1]1 simulates $x \sim \mathcal{S}.q$ and computes an unbiased estimate $\hat{Z}$ of $\frac{\tilde{\pi}(x)}{\mathcal{S}.q(x)}$:

Theorem 1. Given a strategy $\mathcal{S}$ targeting $\pi(x)$, Algorithm [alg:alg1]1 generates $x \sim \mathcal{S}.q$ and a weight $\hat{Z}$ such that $\mathbb{E}[\hat{Z} \mid x] = \frac{\tilde{\pi}(x)}{\mathcal{S}.q(x)} = Z \frac{\pi(x)}{\mathcal{S}.q(x)}$.

When $\mathcal{S}.q$ has a tractable marginal density, Algorithm [alg:alg1]1 computes an exact importance weight. Otherwise, it calls Algorithm [alg:alg2]2, which uses the meta-inference strategy $\mathcal{S}.\mathcal{M}(x)$ to estimate $\frac{1}{\mathcal{S}.q(x)}$. The proof of Theorem 1 is by induction on the level of nesting in the strategy (see Appendix <ref>).

RAVI for MCMC. When models or proposals (or both) in a Metropolis-Hastings sampler do not have tractable closed-form densities, RAVI inference strategies enable computation of MH acceptance probabilities.

We present the algorithm and its justification in Appendix <ref>.

RAVI for Variational Inference. Let $p_\theta(x, y)$ be a latent-variable generative model with parameters $\theta$, and $\mathcal{S}_\theta(y)$ is a family of strategies targeting $p_\theta(x \mid y)$. Given a dataset $y$, variational inference can be applied to maximize (a lower bound on) $\log  p_\theta(y)$, and also to optimize parameters of the posterior approximations in $\mathcal{S}_\theta$, to bring them closer (in KL divergence) to their targets.

\begin{align*}
\LL(p, y, \mathcal{S}) &:= \E [ \log \hat{Z}(p(\cdot, y), \mathcal{S}) ] \le \log p(y)\\ \text{and } \UU(p, y, \mathcal{S}) &:= \E [ \log \check{Z}(p(\cdot, y), \mathcal{S}) ] \ge \log p(y),
\end{align*}

where $\hat{Z}(\tilde{\pi}, \mathcal{S})$ is the estimate returned by  (Alg. [alg:alg1]1) on $\mathcal{S}$ and unnormalized target $\tilde{\pi}$,

and $\check{Z}(\tilde{\pi}, \mathcal{S})$ is the inverse of the weight returned from  (Alg. [alg:alg2]2) when run with unnormalized target $\tilde{\pi}$, inference strategy $\mathcal{S}$, and an exact sample $x \sim \pi$. Because $\hat{Z}$ is an unbiased estimate of $p_\theta(y)$, and ${\check{Z}}^{-1}$ is an unbiased estimate of ${p_\theta(y)}^{-1}$, we have by Jensen's inequality that $\LL(p, y, \mathcal{S})$ and $\UU(p, y, \mathcal{S})$ are lower and upper bounds (respectively) on $\log p_\theta(y)$. As such, we can fit the model parameters $\theta$ to data $y$ by minimizing $\UU(p, y, \mathcal{S})$ or maximizing $\LL(p, y, \mathcal{S})$.

Recursive stochastic gradient estimation.
$\texttt{ELBO}\nabla$ (Alg. [alg:alg3]3) is a procedure for estimating $\LL(p, y, \mathcal{S})$ and its gradient $\grad \LL(p, y, \mathcal{S})$ with respect to the parameters $\theta$ of the model and the strategy. When $(x, y) \sim p(x, y)$, $\texttt{EUBO}\nabla$ (Alg. [alg:alg4]4) estimates $\UU(p, y, \mathcal{S})$ and the gradient $\grad\E_{y \sim p}[\UU(p, y, \mathcal{S})]$.

These procedures employ score function estimation of gradients, but
it is straightforward to incorporate baselines within each procedure to reduce variance.
Depending on $\mathcal{S}$, the reparametrization trick may also be applicable (Appendix <ref>).

Theorem 2.
Given a model $p_\theta(x, y)$ and an inference strategy $\mathcal{S}_\theta$ targeting $p_\theta(x \mid y)$,

Alg. 3 yields unbiased estimates of $\LL(p, y, \mathcal{S})$ and of $\grad \LL(p, y, \mathcal{S})$.
Furthermore, when $(x, y) \sim p_\theta$, Alg. 4 yields
(i) $\hat{U}$ such that
$\E [ \hat{U} \mid y ] = \UU(p, y, \mathcal{S})$,
(ii) $\widehat{\grad}$ such that
$\E [ \widehat{\grad} ] = \grad \E_{y \sim p}[\UU(p, y, \mathcal{S})]$,
and (iii) a value $\mathbf{g}$ such that for any function $R$ that does not depend on $\theta$,
$\E [ \mathbf{g} \cdot R(y) ] = \grad \E_{y \sim p_\theta} [ R(y) ]$  if
$\grad \E_{y \sim p_\theta} [ R(y) ]$ is defined.

In Section 4, we show the tightness of the variational bounds $\LL$ and $\UU$ is given by sums of KL divergences between posterior approximations in $\mathcal{S}_\theta$ and their targets. Thus, optimizing these bounds improves the posterior approximations, either encouraging mass-capturing or mode-seeking behavior.

§ THEORETICAL ANALYSIS

We now present theorems characterizing the quality of RAVI inference: Thm. 3 concerns the variance of weights in a Monte Carlo sampler, and Thm. 4 the tightness of variational bounds. In both cases, error is related to each approximation in the RAVI strategy's divergence to its target posterior.

Sampler variance in Monte Carlo.
Let $\tilde{\pi} = Z\pi$ be an unnormalized target density, and $\mathcal{S}$ an inference strategy targeting $\pi$. As in Section <ref>, we write $\hat{Z}(\tilde{\pi}, \mathcal{S})$ for the weight returned by , and $\text{Var}_{\hat{Z}}(\pi, \mathcal{S})$ for the relative variance of the estimator, $\text{Var}(\frac{\hat{Z}(\tilde{\pi}, \mathcal{S})}{Z})$, which does not depend on $Z$ (and therefore is a function of $\pi$, not $\tilde{\pi}$). Similarly, we write $\check{Z}(\tilde{\pi}, \mathcal{S})$ for the reciprocal of the weight returned by , run with an input $x \sim \pi$. $\text{Var}_{\check{Z}}(\pi, \mathcal{S})$ is its relative variance, $\text{Var}(\frac{Z}{\check{Z}(\tilde{\pi}, \mathcal{S})})$.

Theorem 3. Consider an unnormalized target distribution $\tilde{\pi}(x) = Z\pi(x)$ and an inference strategy $\mathcal{S}$ targeting $\pi(x)$. Then the relative variances of the estimators $\hat{Z}(\tilde{\pi}, \mathcal{S})$ and $\check{Z}(\tilde{\pi}, \mathcal{S})$ are given by the following recursive equations:

\begin{align*}
\text{Var}_{\hat{Z}}&(\pi, \mathcal{S}) = \chi^2(\pi || \mathcal{S}.q) \, +\\
& \mathbb{E}_{x \sim \mathcal{S}.q}\left[\left(\frac{\pi(x)^2}{\mathcal{S}.q(x)^2}\right) \cdot \text{Var}_{\check{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right]\\
\text{Var}_{\check{Z}}&(\pi, \mathcal{S}) = \chi^2(\mathcal{S}.q || \pi) + \\
& \mathbb{E}_{x \sim \pi}\left[\left(\frac{\mathcal{S}.q(x)^2}{\pi(x)^2}\right) \cdot \text{Var}_{\hat{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right]
%\text{Var}\left(\frac{\hat{Z}}{Z}\right) = \sum_{i=0}^{n-1} \mathbb{E}%_{x_0, \dots, x_{i-1}}
%\frac{\hat{Z}_i^2}{Z^2} \cdot \chi^2(&p_{2\lceil \frac{i}{2}\rceil}(x_{i+1} \mid x_{i}, \dots, x_0) \,\, ||\\
%&p_{2\lfloor\frac{i}{2}\rfloor + 1}(x_{i+1} \mid x_i, \dots, x_0))
%w_i^2 \cdot d_i%
\end{align*}

When $\mathcal{S}.q$ is tractable, the second term of each sum is 0.

Tightness of variational bounds.
In VI, the tightness of the variational bounds $\mathcal{L}$ and $\mathcal{U}$ can be characterized as a sum of a KL divergence and a term measuring meta-inference error.

The random variables $\hat{\mathcal{L}}$ and $\hat{\mathcal{U}}$ returned by $\texttt{ELBO}\grad$ and $\texttt{EUBO}\grad$, respectively, are unbiased estimators of $\mathcal{L}(p, y, \mathcal{S})$ and $\mathcal{U}(p, y, \mathcal{S})$, and so can also be viewed as
biased estimators of $\log p(y)$. Writing their bias as $\text{Bias}_\mathcal{L}(p, y, \mathcal{S})$ (and similarly for $\mathcal{U}$), we have:

Theorem 4.
Consider a joint distribution $p(x, y)$ and an inference strategy $\mathcal{S}$ targeting $p(x \mid y)$. Then the following equations give the bias of ${\hat{\LL}}$ and ${\hat{\UU}}$ as estimators of $\log p(y)$:

\begin{align*}
\text{Bias}_\mathcal{L}(p, y, \mathcal{S}) =&\, -\text{KL}(\mathcal{S}.q || p(\cdot \mid y)) \\
&-\mathbb{E}_{x \sim \mathcal{S}.q}[\text{Bias}_\mathcal{U}(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
\text{Bias}_\mathcal{U}(p, y, \mathcal{S}) =&\, \text{KL}(p(\cdot \mid y) || \mathcal{S}.q)\\
&\,\,\,-\mathbb{E}_{x \sim p(\cdot \mid y)}[\text{Bias}_\mathcal{L}(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]
\end{align*}

where the second term in each equation is 0 when $\mathcal{S}.q$ has a tractable marginal density.

Maximizing $\LL$, or minimizing $\UU$, also minimizes these KL divergences. In particular, maximizing $\LL(p, y, \mathcal{S})$ minimizes a `mode-seeking' KL from $\mathcal{S}.q$ to the posterior, whereas minimizing $\mathbb{E}_{y \sim p}[\UU(p, y, \mathcal{S})]$, e.g. by following the gradients of Alg. [alg:alg4]4, implements amortized variational inference, and encourages $\mathcal{S}.q$ to cover the mass of the posterior.

Inference and Meta-Inference. In both Theorems 3 and 4, the first term of the sum is a divergence between $\mathcal{S}.q(x)$, the intractable posterior approximation, and the actual target posterior $p(x \mid y)$. The other term measures the expected quality of meta-inference. Thus the overall error of a RAVI algorithm can be understood as decomposing cleanly into (1) the mismatch between the posterior and the intractable proposal, and (2) the error introduced by meta-inference.

§ EXPERIMENTS

§.§ Improving MCVI

In Section <ref>, we developed a
variant of Salimans et al., 2015's MCVI algorithm
that used SMC for meta-inference.
In Figure <ref>,
we compare vanilla MCVI to the RAVI variant,
with varying $K$ (number of particles used for
meta-inference) and $M$ (number of MCMC steps
in the variational family).

Experimental details. For the MCMC kernel $T$, we use Langevin
ascent with step size $0.015$. For the
meta-inference proposals $R_i(x_{i+1} \rightarrow x_i)$,
we use $\mathcal{N}(x_{i}; f_\mu(x_{i+1}, i), e^{f_{\log \sigma}(x_{i+1}, i)})$, where $f$ is a 4-layer MLP, the step number $i$
is encoded as a one-hot vector, and $f$ outputs the mean $\mu$ and log standard deviation $\log\sigma$ for a conditionally Gaussian proposal.
The same $f$ is used for each experiment, and is trained on
forward rollouts of MCMC (equivalent to using Alg. 3 on $\texttt{rmcvi}$ with $K=1$). The unimodal model is Gaussian with
$\sigma=0.2$, and the multimodal model is a mixture
of 3 Gaussians with standard deviations $0.2, 0.3, \& 2.0$.
The distributions $q_i$ used for importance weighting during sequential Monte Carlo meta-inference are Gaussians with learned $\mu$ and $\sigma$.

Results. Figure <ref> plots the gap $\log p(y) - \mathcal{L}$ for each algorithm's variational bound $\mathcal{L}$. By Theorem 4 this gap is the sum of two terms: $KL(\mathcal{S}.q || p(x \mid y))$ and the expected meta-inference divergence $\mathbb{E}_{x \sim \mathcal{S}.q}[KL(\mathcal{S}.q(x_{0:M} \mid x) || \mathcal{S}.\mathcal{M}(x).q(x_{0:M}))]$. The first term is constant across the algorithms, since they all use the same MCMC-based posterior approximation, so the plots primarily illustrate differences in the quality of meta-inference. MCVI's meta-inference steadily worsens as the chain's length grows, and after 15-25 steps, the meta-inference cost of adding new steps outweighs the benefits to $\mathcal{S}.q$, causing the bound $\mathcal{L}$ to loosen. Our variant, with SMC meta-inference, does not suffer the same penalty, and continues to improve as more steps are added. As discussed in Section <ref>, the same inference strategy ([infstrat:mcvi]) can be used within an importance sampler to derive unbiased marginal likelihood estimates from MCMC runs. The right-hand plot in Figure <ref> shows that this technique can yield accurate estimates with less computation than AIS [Neal, 2001], at least on simple targets. (To fairly account for the computational cost of meta-inference, in the RAVI algorithm we multiply $M$ by $K$ when plotting the total number of MCMC steps.)

Because the variance of AIS is bounded below by sums of divergences between subsequent pairs of intermediate target distributions, the MCMC chain must be long enough to support a very fine annealing schedule, without large jumps. By contrast, RAVI-MCVI requires
only that the marginal distribution of the chain be a good approximation to the posterior, and that SMC meta-inference is sufficiently accurate. For some problems, this may be less expensive than the long chain required by AIS.

Improving Markov Chain Variational Inference with RAVI. Left and Middle: On unimodal and multimodal targets, MCVI begins to degrade after 15-25 steps of MCMC. RAVI-MCVI with sufficiently many particles continues to improve as more MCMC steps are added.

Right: When MCMC converges quickly to a reasonable approximation of the posterior, RAVI-MCVI can give more accurate estimates of marginal likelihoods than standard techniques such as AIS. The $x$ axis of this plot counts total MCMC steps simulated, whether as part of inference or meta-inference; for RAVI-MCVI($K$), this is $KM$, where $M$ is the length of the forward Markov chain and $K$ is the number of SMC particles used for meta-inference.

§.§ Agglomerative Clustering for Dirichlet Process Mixtures

A promising application of RAVI is to transform heuristic randomized algorithms into
unbiased and consistent Monte Carlo estimators, by using them as proposal distributions. In this section, we design a RAVI inference strategy for clustering in Dirichlet process mixtures, based on a randomized agglomerative clustering algorithm (Inference Strategy [infstrat:agglom]2).

RAVI Inference Strategy 2: Agglomerative Clustering
agglomagglom($X, K$).qagglommagglom($X, K$).M($\Pi$).q
infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

partition $\Pi$ of dataset $X$
merge sequence $\textit{mv}_{1:|X|-|\Pi|}$
$\Pi \gets \{\{x\} \mid x \in X\}$ *Initial partition
$l \in 1, \dots, |X|$
unordered pair $\{i, j\}$ of clusters in $\Pi$
$w_{\{i,j\}} \gets \pi((\Pi \setminus \{i, j\}) \cup \{i \cup j\})$

$w_{\text{stop}} = \pi(\Pi)$
$\textit{All} \gets \{\text{stop}\} \cup \{\{i, j\} \mid i, j \in \Pi\}$
$\textit{mv}_l \sim \text{Discrete}(\{m \Rightarrow w_m \mid m \in \textit{All}\})$
$\textit{mv}_l = \text{stop}$$\textbf{break}$
$\Pi \gets (\Pi \setminus \textit{mv}_l) \cup (\cup \textit{mv}_l)$*Perform the merge

metaalgMeta-Posterior Approx.

merge sequence $\textit{mv}_{1:|X|-|\Pi|}$
particles $\textit{mv}_{1:|X|-|\Pi|}^{1:K}$, ancestors $a_{1:|X|-|\Pi|}^{1:K}$
$k \in 1,\dots, K$$\Pi_0^k, \textit{tr}_k  \gets \{\{x\} \mid x \in X\}, []$
$l \in 1, \dots, |X|-|\Pi|$
$k \in 1,\dots,K$
unordered pair $\{i, j\}$ in $\Pi_{l-1}^k$
$w_{\{i,j\}} \gets \pi((\Pi_{l-1}^k \setminus \{i, j\}) \cup \{i \cup j\})$

$w_{\text{stop}} = \pi(\Pi_{l-1}^k)$
$\textit{All} \gets \{\text{stop}\} \cup \{\{i, j\} \mid i, j \in \Pi_{l-1}^k\}$
$\textit{Ok} \gets \{\{i, j\} \in \textit{All} \mid \exists c \in \Pi. i \cup j \subseteq c\}$
$\textit{mv}_l \sim \text{Discrete}(\{m \Rightarrow w_m \mid m \in \textit{Ok}\})$
$\Pi_{l-1}^k \gets (\Pi_{l-1}^k \setminus \textit{mv}_l) \cup (\cup \textit{mv}_l)$
$\textit{tr}_k \gets [\textit{tr}_k\dots, \textit{mv}_l]$
$W_l^k \gets \frac{\sum_{m \in \textit{Ok}} w_m}{\sum_{m \in \textit{All}} w_m}$

$k \in 1,\dots,K$
$a_l^k \sim \text{Discrete}(W_l^{1:K})$
*resampling step
$\Pi_l^k, \textit{tr}_k \gets \Pi_{l-1}^{a_l^k}, \textit{tr}_{a_l^k}$

$[\textit{tr}_1\dots, \text{stop}]$
agglommmagglom($X, K$).M($\Pi$).M($\textit{mv}_{1:|X|-|\Pi|}$).q

Conditional SMC (omitted for space, but similar to that of )

Datasets and Models. We test our algorithm on three clustering problems. The first is a synthetic 1D dataset sampled from a Dirichlet process (DP) mixture prior. The second is a standard benchmark dataset of galaxy velocities [Drinkwater et al., 2004, Cusumano-Towner and Mansinghka, 2017], which we model using a DP mixture with Gaussian likelihoods and $\alpha = 1$.  The last is a data-cleaning task, correcting typos in 1k strings from Medicare records [Medicare, 2012]. We adapt the generative model of Lew et al., 2021. Using an English character-level bigram model $H(s) = h(s_1) \prod_{i=2}^{|s|} h(s_i \mid s_{i-1})$, we model the data $\{y_i\}$  with a DP prior:

\begin{align*}
G &\sim DP(H, \alpha = 1.0)\\
x_i \mid G &\sim G\\
y_i \mid x_i &\sim f(\cdot \mid x_i)
\end{align*}

Here, the likelihood $f(y_i \mid x_i)$ models typos. We set $f$ to be

\[
f(y_i \mid x_i) \propto
\begin{cases}
\mathbf{1}[x_i = y_i] & (x_i, y_i) \not\in \mathcal{L} \times \mathcal{L} \\
\frac{\text{NegBin}(\tau(x_i, y_i); \lceil \frac{|s|}{5} \rceil, 0.9)}{(5.09|s|)^{\tau(x_i, y_i)}} & (x_i, y_i) \in \mathcal{L} \times \mathcal{L}
\end{cases},
\]

where $\tau(x_i, y_i)$ is the Damerau-Levenshtein edit distance between $x_i$ and $y_i$, and $\mathcal{L}$ is the set of all observed strings $\{y \mid \exists i. \, y = y_i\}$.
[We assume that the data $\mathcal{L}$ includes at least one example of every clean string.

When $x_i \in \mathcal{L}$, we model a negative-binomially distributed number of typos, where the number of trials depends on the length of the string.]
We perform inference in a collapsed version of the model, with the $x_i$ marginalized out:

\begin{align*}
\Pi &\sim CRP(n = N, \alpha = 1.0)\\
y_I \mid \Pi &\sim F(y_I).
\end{align*}

Here $\Pi$ is a partition, $I$ ranges over the components of $\Pi$ (each of which is a subset of data indices), and $F(y_I) = \sum_{x \in \mathcal{L}} h(x) \prod_{i \in I} f(y_i \mid x)$ is the marginal likelihood of $y_I$ as a sequence of noisy observations of a latent string.

Baseline. We compare to an SMC baseline, inspired by PClean's inference [Lew et al., 2021], that targets a sequence of posteriors, where the $t^{\text{th}}$ posterior incorporates the first $t$ datapoints.
The SMC proposal is locally optimal, assigning the newest datapoint to an existing component $I$ with probability proportional to $\frac{|I|}{t + \alpha - 1} \cdot F(y_{I} \cup \{y_t\})$, or to a new component with probability proportional to $\frac{\alpha}{t + \alpha - 1} \cdot F(\{y_t\})$. We do not compare to a Gibbs sampling baseline, as Gibbs sampling does not yield marginal likelihood estimates, but do perform a Gibbs rejuvenation sweep every 20 iterations of SMC.

RAVI algorithm. We apply Algorithm [alg:alg1]1 to the inference strategy $\texttt{agglom}$ (Inference Strategy [infstrat:agglom]2). The strategy is based
on a randomized agglomerative clustering algorithm: each datapoint begins
in its own cluster (L1), and we repeatedly
choose to either merge two clusters (L9) or stop and propose the current partition (L8). The sequence of merge
decisions $\textit{mv}_1, \dots, \textit{mv}_{|X|-|\Pi|}$ are the auxiliary variables
of our proposal distribution; the final output is the proposed clustering $\Pi$. Our meta-inference $\texttt{agglom}(X, K).\mathcal{M}(\Pi).q$ infers the sequence of merges from the observed clustering $\Pi$, using $K$-particle SMC with proposals that mimic the forward process but choose only from a restricted set $\textit{Ok}$ of possible merges (L8), to avoid making any choices that disagree with $\Pi$. SMC introduces additional auxiliary variables, so we also include a conditional SMC meta-meta-posterior approximation (not shown, but nearly identical to [infstrat:mcvi]'s).

Results. Table <ref> shows
average log marginal likelihood estimates; higher is
better. On synthetic Gaussian data, the algorithms
perform comparably. On the galaxy data, RAVI
agglomerative clustering finds modes that SMC misses,
leading to a 3-nat improvement in the
average log marginal likelihood. In the Medicare
data example, SMC misses the ground-truth clustering
and hypothesizes many unlikely typos to
explain the data. The RAVI agglomerative clustering is
less greedy, considering $O(N^2)$ possible merges at each
step, rather than $O(N)$. As such, it is able to
find the ground truth clustering, correctly identifying
all typos (unlike PClean [Lew et al., 2021], which achieves
only 90% accuracy on this dataset) and reporting a log marginal likelihood thousands of
nats higher than the SMC algorithm.

2l|Gaussian likelihood [Cusumano-Towner and Mansinghka, 2017], synthetic data

SMC + adapted proposals        $-125.09 \pm 0.38 $

RAVI agglomerative clustering                        ${-125.97 \pm 1.62}$

2l|Gaussian likelihood [Cusumano-Towner and Mansinghka, 2017], Galaxy data [Drinkwater et al., 2004]

SMC + adapted proposals                      $-426.20 \pm 1.26  $

RAVI agglomerative clustering                        $\mathbf{-423.03 \pm 0.94}$

2l|PClean typos likelihood [Lew et al., 2021], Hospital data [Medicare, 2012]

SMC + adpated proposals                $-40,239 \pm 1,532$

RAVI agglomerative clustering                          $\mathbf{-13,851.0 \pm 0.01}  $

RAVI agglomerative clustering vs. SMC baseline.

§ RELATED WORK AND DISCUSSION

Related work. RAVI builds on and generalizes
recent work from both the Monte Carlo and variational
inference literatures. For example, Salimans et al., 2015 and Ranganath et al., 2016 showed how
auxiliary variables could be used to
construct and optimize variational bounds for
specific families of expressive variational approximations.
Sobolev and Vetrov, 2019 presented tighter bounds
in a more general setting. RAVI is a further
generalization, in two directions: first, we show
that these bounds arise from particular
choices of meta-inference strategy, and can be
tightened by improving meta-inference; and second, we
extend the results to the Monte Carlo setting, enabling
learned variational families to be used as IS, SMC, or MH proposals. We also provide general theorems about the
variance of RAVI samplers and the bias of RAVI variational bounds, which can be applied to analyze both new and existing algorithms.

RAVI is also related to other compositional or unifying frameworks for thinking about broad classes of inference algorithms [Mansinghka et al., 2014, Zinkov and Shan, 2016, Ścibior et al., 2017, Ścibior et al., 2018, Mansinghka et al., 2018, Cusumano-Towner et al., 2019, Stites et al., 2021, Neklyudov et al., 2020, Zimmermann et al., 2021, Andrieu et al., 2020, Storvik, 2011, Finke, 2015, Finke and Thiery, 2019], some of which involve recursive constructions [Naesseth et al., 2015, Del Moral et al., 2006, Domke and Sheldon, 2019]. However, to our knowledge, RAVI's inference strategies are novel. For example, although RAVI and Nested IS (NIS) [Naesseth et al., 2015] are both approaches to inference with `intractable proposals,' NIS approximately samples a proposal distribution with a tractable (unnormalized) density, whereas RAVI approximates the density of a proposal that can be simulated tractably, but whose marginal density (even unnormalized) is intractable.
As another example, Domke and Sheldon, 2019's framework of estimator-coupling pairs constructs variational bounds and marginal likelihood estimators recursively, but unlike in RAVI,
the posterior approximations cannot be used to formulate objectives for amortized VI or as components of Metropolis-Hastings proposals.

Finally, researchers have used meta-inference to construct bounds on KL divergences [Cusumano-Towner and Mansinghka, 2017] and other information-theoretic quantities [Saad et al., 2022]. In Appendix <ref>, we show how to apply such bounds in the general RAVI setting.

Outlook and Limitations. RAVI expands the design space for Monte Carlo and variational inference. It gives unifying correctness proofs for over a dozen methods from the literature, and novel theorems that  characterize their behavior. Experiments show that RAVI helps to design algorithms that significantly improve accuracy over previously introduced Monte Carlo and variational inference methods. However, some difficulties remain. For example, the gradient estimators we present (Algs. [alg:alg3]3 and [alg:alg4]4)
have high variance for some strategies $\mathcal{S}$; in Appendix <ref>, we give estimators that exploit the reparameterization trick, but they only help when the proposals in $\mathcal{S}$ can be
reparameterized, which is not the case, e.g., for SMC. In these cases, RAVI can still be used to derive
objectives for optimization, but practitioners will need other ways of reducing the variance of gradient estimates; many results from
the literature [Mnih and Rezende, 2016, Tucker et al., 2018] should apply.

Another difficulty is that RAVI algorithms can be complex to implement. We are exploring an automated implementation based on PPLs [Cusumano-Towner et al., 2019, van de Meent et al., 2018]: if the posterior and meta-posterior approximations in a RAVI strategy $\mathcal{S}$ are probabilistic programs, we can provide Algs. 1-4 as higher-order functions, which automate the necessary densities and gradients.

[Andrieu et al., 2010]
Christophe Andrieu, Arnaud Doucet, and Roman Holenstein.
Particle markov chain monte carlo methods.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 720 (3):0 269–342, 2010.

[Andrieu et al., 2018]
Christophe Andrieu, Arnaud Doucet, Sinan Yıldırım, and Nicolas Chopin.
On the utility of metropolis-hastings with asymmetric acceptance
arXiv preprint arXiv:1803.09527, 2018.

[Andrieu et al., 2020]
Christophe Andrieu, Anthony Lee, and Sam Livingstone.
A general perspective on the metropolis-hastings kernel.
arXiv preprint arXiv:2012.14881, 2020.

[Bachman and Precup, 2015]
Philip Bachman and Doina Precup.
Training deep generative models: Variations on a theme.
In NIPS Approximate Inference Workshop, 2015.

[Burda et al., 2015]
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov.
Importance weighted autoencoders.
arXiv preprint arXiv:1509.00519, 2015.

[Chatterjee et al., 2018]
Sourav Chatterjee, Persi Diaconis, et al.
The sample size required in importance sampling.
The Annals of Applied Probability, 280 (2):0
1099–1135, 2018.

[Chopin and Papaspiliopoulos, 2020]
Nicolas Chopin and Omiros Papaspiliopoulos.
An introduction to sequential Monte Carlo.
Springer, 2020.

[Chopin et al., 2013]
Nicolas Chopin, Pierre E Jacob, and Omiros Papaspiliopoulos.
Smc2: an efficient algorithm for sequential analysis of state space
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 750 (3):0 397–426, 2013.

[Cremer et al., 2017]
Chris Cremer, Quaid Morris, and David Duvenaud.
Reinterpreting importance-weighted autoencoders.
arXiv preprint arXiv:1704.02916, 2017.

[Cusumano-Towner and Mansinghka, 2017]
Marco F Cusumano-Towner and Vikash K Mansinghka.
Aide: An algorithm for measuring the accuracy of probabilistic
inference algorithms.
arXiv preprint arXiv:1705.07224, 2017.

[Cusumano-Towner et al., 2019]
Marco F Cusumano-Towner, Feras A Saad, Alexander K Lew, and Vikash K
Gen: a general-purpose probabilistic programming system with
programmable inference.
In Proceedings of the 40th ACM SIGPLAN Conference on
Programming Language Design and Implementation, pages 221–236, 2019.

[Del Moral et al., 2006]
Pierre Del Moral, Arnaud Doucet, and Ajay Jasra.
Sequential monte carlo samplers.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 680 (3):0 411–436, 2006.

[Domke and Sheldon, 2019]
Justin Domke and Daniel Sheldon.
Divide and couple: Using monte carlo variational objectives for
posterior approximation.
arXiv preprint arXiv:1906.10115, 2019.

[Doucet et al., 2000]
Arnaud Doucet, Simon Godsill, and Christophe Andrieu.
On sequential monte carlo sampling methods for bayesian filtering.
Statistics and computing, 100 (3):0 197–208,

[Drinkwater et al., 2004]
Michael J Drinkwater, Quentin A Parker, Dominique Proust, Eric Slezak, and
Hernán Quintana.
The large scale distribution of galaxies in the shapley supercluster.
Publications of the Astronomical Society of Australia,
210 (1):0 89–96, 2004.

[Finke, 2015]
Axel Finke.
On extended state-space constructions for Monte Carlo methods.
PhD thesis, University of Warwick, 2015.

[Finke and Thiery, 2019]
Axel Finke and Alexandre H Thiery.
On importance-weighted autoencoders.
arXiv preprint arXiv:1907.10477, 2019.

[Lew et al., 2021]
Alexander Lew, Monica Agrawal, David Sontag, and Vikash Mansinghka.
Pclean: Bayesian data cleaning at scale with domain-specific
probabilistic programming.
In International Conference on Artificial Intelligence and
Statistics, pages 1927–1935. PMLR, 2021.

[Liu, 2008]
Jun S Liu.
Monte Carlo strategies in scientific computing.
Springer Science & Business Media, 2008.

[Mansinghka et al., 2014]
Vikash Mansinghka, Daniel Selsam, and Yura Perov.
Venture: a higher-order probabilistic programming platform with
programmable inference.
arXiv preprint arXiv:1404.0099, 2014.

[Mansinghka et al., 2018]
Vikash K Mansinghka, Ulrich Schaechtle, Shivam Handa, Alexey Radul, Yutian
Chen, and Martin Rinard.
Probabilistic programming with programmable inference.
In Proceedings of the 39th ACM SIGPLAN Conference on
Programming Language Design and Implementation, pages 603–616, 2018.

[Medicare, 2012]
Hospital compare, 2012.

[Mnih and Rezende, 2016]
Andriy Mnih and Danilo Rezende.
Variational inference for monte carlo objectives.
In International Conference on Machine Learning, pages
2188–2196. PMLR, 2016.

[Murphy, 2012]
Kevin P Murphy.
Machine learning: a probabilistic perspective.
MIT press, 2012.

[Naderiparizi et al., 2019]
Saeid Naderiparizi, Adam Ścibior, Andreas Munk, Mehrdad Ghadiri,
Atılım Güneş Baydin, Bradley Gram-Hansen, Christian Schroeder de Witt,
Robert Zinkov, Philip H. S. Torr, Tom Rainforth, Yee Whye Teh, and Frank
Amortized rejection sampling in universal probabilistic programming,

[Naesseth et al., 2015]
Christian Naesseth, Fredrik Lindsten, and Thomas Schon.
Nested sequential monte carlo methods.
In International Conference on Machine Learning, pages
1292–1301. PMLR, 2015.

[Naesseth et al., 2018]
Christian Naesseth, Scott Linderman, Rajesh Ranganath, and David Blei.
Variational sequential monte carlo.
In International Conference on Artificial Intelligence and
Statistics, pages 968–977. PMLR, 2018.

[Neal, 2001]
Radford M Neal.
Annealed importance sampling.
Statistics and computing, 110 (2):0 125–139,

[Neklyudov et al., 2020]
Kirill Neklyudov, Max Welling, Evgenii Egorov, and Dmitry Vetrov.
Involutive mcmc: a unifying framework.
In International Conference on Machine Learning, pages
7273–7282. PMLR, 2020.

[Newton and Raftery, 1994]
Michael A Newton and Adrian E Raftery.
Approximate bayesian inference with the weighted likelihood
Journal of the Royal Statistical Society: Series B
(Methodological), 560 (1):0 3–26, 1994.

[Ranganath et al., 2016]
Rajesh Ranganath, Dustin Tran, and David Blei.
Hierarchical variational models.
In International conference on machine learning, pages
324–333. PMLR, 2016.

[Saad et al., 2022]
Feras A Saad, Marco Cusumano-Towner, and Vikash K Mansinghka.
Estimators of entropy and information via inference in probabilistic
arXiv preprint arXiv:2202.12363, 2022.

[Salimans et al., 2015]
Tim Salimans, Diederik Kingma, and Max Welling.
Markov chain monte carlo and variational inference: Bridging the gap.
In International Conference on Machine Learning, pages
1218–1226. PMLR, 2015.

[Ścibior et al., 2017]
Adam Ścibior, Ohad Kammar, Matthijs Vákár, Sam Staton, Hongseok
Yang, Yufei Cai, Klaus Ostermann, Sean K Moss, Chris Heunen, and Zoubin
Denotational validation of higher-order bayesian inference.
arXiv preprint arXiv:1711.03219, 2017.

[Ścibior et al., 2018]
Adam Ścibior, Ohad Kammar, and Zoubin Ghahramani.
Functional programming for modular bayesian inference.
Proceedings of the ACM on Programming Languages, 20
(ICFP):0 1–29, 2018.

[Sobolev and Vetrov, 2019]
Artem Sobolev and Dmitry P Vetrov.
Importance weighted hierarchical variational inference.
Advances in Neural Information Processing Systems, 32, 2019.

[Stites et al., 2021]
Sam Stites, Heiko Zimmermann, Hao Wu, Eli Sennesh, et al.
Learning proposals for probabilistic programs with inference
arXiv preprint arXiv:2103.00668, 2021.

[Storvik, 2011]
Geir Storvik.
On the flexibility of metropolis–hastings acceptance probabilities
in auxiliary variable proposal generation.
Scandinavian Journal of Statistics, 380 (2):0
342–358, 2011.

[Tucker et al., 2018]
George Tucker, Dieterich Lawson, Shixiang Gu, and Chris J Maddison.
Doubly reparameterized gradient estimators for monte carlo
arXiv preprint arXiv:1810.04152, 2018.

[van de Meent et al., 2018]
Jan-Willem van de Meent, Brooks Paige, Hongseok Yang, and Frank Wood.
An introduction to probabilistic programming.
arXiv preprint arXiv:1809.10756, 2018.

[Zimmermann et al., 2021]
Heiko Zimmermann, Hao Wu, Babak Esmaeili, and Jan-Willem van de Meent.
Nested variational inference.
arXiv preprint arXiv:2106.11302, 2021.

[Zinkov and Shan, 2016]
Robert Zinkov and Chung-chieh Shan.
Composing inference algorithms as program transformations.
arXiv preprint arXiv:1603.01882, 2016.

§ SUPPLEMENTARY MATERIAL FOR “RECURSIVE MONTE CARLO AND VARIATIONAL INFERENCE WITH AUXILIARY VARIABLES”

This document and the accompanying code files contain supplementary material for the submission “Recursive Monte Carlo and Variational Inference with Auxiliary Variables.” In particular, we provide:

* In Section <ref>, proofs of Theorems 1-4.

* In Section <ref>, RAVI inference strategies for many existing algorithms.

* In Section <ref>, a further discussion of the absolute continuity requirements for RAVI and how they can be relaxed.

* In Section <ref>, other applications of RAVI inference strategies, to parameterize rejection sampling and KL divergence estimation algorithms.

§ OMITTED PROOFS.

Throughout this section, we use the notation introduced in Section <ref>: the random variable $\hat{Z}(\tilde{\pi}, \mathcal{S})$ is the weight returned by $\texttt{IMPORTANCE}(\tilde{\pi}, \mathcal{S})$, and $\check{Z}(\tilde{\pi}, \mathcal{S})$ is the reciprocal of the weight returned by $\texttt{HME}(\tilde{\pi}, x, \mathcal{S})$, for $x \sim \pi$.

§.§ Proof of Theorem 1.

We prove a stronger version of Theorem 1 (in order to use a stronger inductive hypothesis in the proof):

Theorem 1.
Let $\tilde{\pi}(x) = Z\pi(x)$ be an unnormalized target density, and $\mathcal{S}$ an inference strategy targeting $\pi(x)$. Then:

* $\texttt{IMPORTANCE}(\mathcal{S}, \tilde{\pi})$ generates $(x, \hat{Z})$ with $x \sim \mathcal{S}.q$ and $\mathbb{E}[\hat{Z} \mid x] = Z\frac{\pi(x)}{\mathcal{S}.q(x)}$. Furthermore, the unconditional expectation $\mathbb{E}[\hat{Z}(\tilde{\pi}, \mathcal{S})] = Z$.

* $\mathbb{E}[{\check{Z}}(\tilde{\pi}, \mathcal{S})^{-1}] = \mathbb{E}_{x \sim \pi}[\texttt{HME}(\mathcal{S}, x, \tilde{\pi})] = Z^{-1}.$

Proof. The proof is by induction on the level of nesting
present in the inference strategy.

First consider the case where $\mathcal{S}.q$ has a tractable
marginal density. Then:

*  samples $x \sim \mathcal{S}.q$ on line 2, and
computes $\hat{Z} = \frac{\tilde{\pi}(x)}{\mathcal{S}.q(x)} = Z\frac{\pi(x)}{\mathcal{S}.q(x)}$ exactly (lines 3 and 7). By the standard importance sampling argument, the unconditional expectation $\mathbb{E}[\hat{Z}(\tilde{\pi}, \mathcal{S})] = \mathbb{E}_{x \sim \mathcal{S}.q}[Z\frac{\pi(x)}{\mathcal{S}.q(x)}]=Z\mathbb{E}_{x \sim \pi}[1] = Z$. (This argument relies on the fact that, because $\mathcal{S}$ targets $\pi$, $\pi$ is absolutely continuous with respect to $\mathcal{S}.q$.)

* $\texttt{HME}(\mathcal{S}, x, \tilde{\pi})$ returns exactly $\frac{\mathcal{S}.q(x)}{\tilde{\pi}(x)}$ (lines 2 and 5), and
$$\mathbb{E}_{x \sim \pi}\left[\frac{\mathcal{S}.q(x)}{\tilde{\pi}(x)}\right] = \int \pi(x) \frac{\mathcal{S}.q(x)}{Z\pi(x)}\text{d}x = \frac{1}{Z} \int \mathcal{S}.q(x) \text{d}x = \frac{1}{Z},$$
where the last step follows because $\mathcal{S}.q$ is a normalized probability density, and $\mathcal{S}.q$ is absolutely continuous with respect to $\pi$.

Now consider the inductive step. Assume $\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x)\text{d}r$ and that for all $x$, the theorem holds for the inference strategy $\mathcal{S}.\mathcal{M}(x)$ and the unnormalized target distribution $\mathcal{S}.q(\cdot, x)$. In this case:

* On line 5,  generates $x \sim \mathcal{S}.q$ and $r \sim \mathcal{S}.q(r \mid x)$. In the call to $\texttt{HME}$, the unnormalized target distribution is $\mathcal{S}.q(\cdot, x)$, and so the normalizing constant is $\mathcal{S}.q(x)$ and the normalized target is $\mathcal{S}.q(r \mid x)$. By the inductive hypothesis, the call to $\texttt{HME}$ on line 6 returns an unbiased estimate of the normalizing constant's reciprocal, i.e. $\mathbb{E}[w \mid x] = \frac{1}{\mathcal{S}.q(x)}$. Since  returns $\hat{Z} = w\tilde{\pi}(x)$ on line 7, this implies that $\mathbb{E}[\hat{Z} \mid x] = \frac{\tilde{\pi}(x)}{\mathcal{S}.q(x)} = Z\frac{\pi(x)}{\mathcal{S}.q(x)}$. From this, the same standard importance sampling argument as above shows that the unconditional expectation $\mathbb{E}[\hat{Z}(\tilde{\pi}, \mathcal{S})] = Z$.

* On line 4,  calls  on the
unnormalized target $\mathcal{S}.q(\cdot, x)$, and so by the inductive hypothesis, $\mathbb{E}[w] = \mathcal{S}.q(x)$ (the normalizing constant). On line 5, the returned weight has expectation $\mathbb{E}_{x \sim \pi}\left[\frac{w}{\tilde{\pi}(x)}\right]=\frac{1}{Z} \int \pi(x) \cdot \frac{\mathcal{S}.q(x)}{\pi(x)} \text{d}x = \frac{1}{Z}$, where the last equality again follows because $\mathcal{S}.q$ is a normalized density, and $\mathcal{S}.q$ is absolutely continuous with respect to $\pi$.

§.§ Proof of Theorem 2

For an inference strategy $\mathcal{S}$ targeting $p(x \mid y)$, if $\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x)\text{d}r$ has an intractable marginal density, then:

\[
\LL(p, y, \mathcal{S}) = \E_{x \sim \mathcal{S}.q} [\log p(x, y) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)]
\]

\[
\UU(p, y, \mathcal{S}) = \E_{x \sim p( \cdot | y) } [ \log p(x, y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)) ]
\]

Proof. For the first conclusion,

\begin{align}
\LL(p, y, \mathcal{S})
&= \E[\log \hat{Z}(p(\cdot \mid y), \mathcal{S})]\\
&= \E\left[\log \frac{p(x, y)}{\check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))}\right]\\
&= \E_{x \sim \mathcal{S}.q}[\E[\log p(x, y) - \log \check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x)) \mid x]]\\
&= \E_{x \sim \mathcal{S}.q}[\log p(x, y) - \mathbb{E}[\log \check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x)) \mid x]]\\
&= \E_{x \sim \mathcal{S}.q} [\log p(x, y) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)) ]
\end{align}

The same approach, but with $\E [ \log \check{Z} ]$, can be used to prove the other conclusion.

Theorem 2.
Given a model $p_\theta(x, y)$ and an inference strategy $\mathcal{S}_\theta$ targeting $p_\theta(x \mid y)$,

Alg. 3 yields unbiased estimates of $\LL(p, y, \mathcal{S})$ and of $\grad \LL(p, y, \mathcal{S})$.
Furthermore, when $(x, y) \sim p_\theta$, Alg. 4 yields
(i) $\hat{U}$ such that
$\E [ \hat{U} \mid y ] = \UU(p, y, \mathcal{S})$,
(ii) $\widehat{\grad}$ such that
$\E [ \widehat{\grad} ] = \grad \mathbb{E}_{y \sim p_\theta}[\UU(p, y, \mathcal{S})]$,
and (iii) a value $\mathbf{g}$ such that for any function $R$ that does not depend on $\theta$,
$\E [ \mathbf{g} \cdot R(y) ] = \grad \E_{y \sim p_\theta} [ R(y) ]$  if
$\grad \E_{y \sim p_\theta} [ R(y) ]$ is defined.

Proof. The proof is by induction on the level of nesting present in the inference strategy.

First consider inference strategies $\mathcal{S}$ with tractable proposals $\mathcal{S}.q(x)$. In this case $\texttt{ELBO}\nabla$ generates $x \sim \mathcal{S}.q$ and returns $\hat{L} = \log p(x, y) - \log \mathcal{S}.q(x)$ and $\widehat{\grad} = \grad (\log p(x, y) - \log \mathcal{S}.q(x)) + (\grad \log \mathcal{S}.q(x))(\log p(x, y) - \log \mathcal{S}.q(x))$. Clearly, $\E_{x \sim \mathcal{S}.q}[\hat{L}] = \E[\log \hat{Z}(p(\cdot, y), \mathcal{S})] = \LL(p, y, \mathcal{S})$. And by the log-derivative trick, $\E_{x \sim \mathcal{S}.q}[\widehat{\grad}] = \E[\grad (\log p(x, y) - \log \mathcal{S}.q(x))] = \E[\mathcal{L}(p, y, \mathcal{S})]$. When we apply $\texttt{EUBO}\nabla$ to $\mathcal{S}$ with $(x, y) \sim p$, it returns (1) $\hat{U} = \log p(x, y) - \log \mathcal{S}.q(x)$ (for which $\E[\hat{U} \mid y] = \UU(p, y, \mathcal{S})$), (2) $\widehat{\grad} = \grad (\log p(x, y) - \log \mathcal{S}.q(x)) + \grad \log p(x, y) (\log p(x, y) - \log \mathcal{S}.q(x))$ (for which, by the log-derivative trick, $\mathbb{E}[\widehat{\grad}] = \grad \E_{y \sim p}[\UU(p, y, \mathcal{S})]$), and (3) $\mathbf{g} = \grad \log p(x, y)$. This last return value satisfies the spec for $\mathbf{g}$ because if $R$ does not depend on $\theta$, then $\mathbb{E}_{(x, y) \sim p}[R(y) \cdot \grad\log p(x, y)] = \int\int p(x, y) \cdot
\frac{\grad p(x, y)}{p(x, y)} \cdot R(y) \text{d}x \text{d}y = \grad \int \int p(x, y) R(y) \text{d}x \text{d}y = \grad \E [R(y)]$, as required.

Now consider the inductive step. Assume the theorem holds for the inference strategy $\mathcal{S}.\mathcal{M}(x)$ and joint distribution $\mathcal{S}.q(r, x)$.

We first consider $\texttt{ELBO}\nabla$. It generates $(r, x) \sim \mathcal{S}.q$ before calling $\texttt{EUBO}\nabla$, which by induction returns $(\hat{U}, \widehat{\grad}, \mathbf{g})$ such that:

* $\mathbb{E}[\hat{U} \mid x] = \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))$

* $\E[\widehat{\grad}] = \grad\E_{x \sim \mathcal{S}.q}[\UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]$

* $\E[g \cdot R(x)] = \grad \E_{x \sim \mathcal{S}.q}[R(x)]$ for all valid $R$.

$\texttt{ELBO}\nabla$ computes its first return value, $\hat{L}$, as $\log p(x, y) - \hat{U}$, so

\begin{align*}
\E[\hat{L}] &= \E[\log p(x, y) - \hat{U}]\\
&= \E_{x \sim \mathcal{S}.q}[\E[\log p(x, y) - \hat{U} \mid x]]\\
&= \E_{x \sim \mathcal{S}.q}[\log p(x, y) - \E[\hat{U} \mid x]]\\
&= \E_{x \sim \mathcal{S}.q}[\log p(x, y)-\UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= \LL(p, y, \mathcal{S}),
\end{align*}

where the fourth equality holds by the inductive hypothesis and the final one by Lemma 1. Its second return value is computed as $\widehat{\grad}' = \grad\log p(x, y) + \mathbf{g}\log p(x, y) - \widehat{\grad}$, and so

\begin{align*}
\E [ \widehat{\grad}' ]
&= \E \left[ \grad \log p(x, y) + \mathbf{g} \cdot \log p(x, y) - \widehat{\grad} \right]\\
&= \E \left[ \grad \log p(x, y) \right]
+ \grad \E_{x \sim \mathcal{S}.q} [ \log \boxed{p}(x, y) ]
- \grad \E_{x \sim \mathcal{S}.q} [ \mathcal{U}(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)) ]\\
&= \grad \E \left[ \log p(x, y) \right]
- \grad \E_{x \sim \mathcal{S}.q} [ \mathcal{U}(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)) ]\\
&= \grad \E_{x \sim \mathcal{S}.q} [ \log p(x, y) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)) ]\\
&= \grad \LL(p, y, \mathcal{S}),
\end{align*}

where $\boxed{p}(x, y)$ denotes the distribution $p(x, y)$ but without a dependence on $\theta$, for the purposes of differentiation with respect to $\theta$. The second equality holds by the inductive hypothesis about $\mathbf{g}$ (with $R(x) = \log \boxed{p}(x, y)$) and about $\widehat{\grad}$, and the third uses the log-derivative trick. The final equation is due to Lemma 1.

We now turn to $\texttt{EUBO}\nabla$.
By induction, the call to $\texttt{ELBO}\nabla$ satisfies the theorem, and so:

* $\E [ \hat{L} ] = \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))$

* $\E [ \widehat{\grad} ] = \grad\LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))$

We treat each of the return values, $(\hat{U}, \widehat{\grad}, \mathbf{g})$, in sequence. We view them as random variables, accounting for stochasticity in the algorithm as well as the inputs $(x, y)$, which are assumed in the theorem's statement to be jointly distributed according to $p$.

First, $\hat{U}$ is computed as $\log p(x, y) - \hat{L}$, and so

\begin{align*}
\E [ \hat{U} | y ]
&= \E_{x \sim p(\cdot | y)} [ \E [ \log p(x, y) - \hat{L} | x, y ] ]\\
&= \E_{x \sim p(\cdot | y) } [ \log p(x, y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)) ]\\
&= \UU(p, y, \mathcal{S}).
\end{align*}

Next, $\E[ \widehat{\grad}']$:

\begin{align*}
\E[\widehat{\grad}']
&= \E_{x, y \sim p} \left[
\E \left[
\grad \log p(x, y)
+ (\grad \log p(x, y)) \cdot \hat{U}
- \widehat{\grad} | x, y
\right]
\right]\\
&= \E_{x, y \sim p} \left[
\grad \log p(x, y) +
(\grad \log p(x, y))
\cdot \E \left[ \hat{U} | x, y \right]
- \E \left[\widehat{\grad} | x, y \right]
\right]\\
&= \E_{x, y \sim p} \left[
\grad \log p(x, y) +
(\grad \log p(x, y))
\cdot \E \left[ \log p(x, y) - \hat{L} \mid x, y \right]
- \grad \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))
\right]\\
&= \E_{x, y \sim p} \left[
\grad \log p(x, y)
- \grad \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))
+ (\grad \log p(x, y))
\cdot (\log p(x, y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))
\right]\\
&= \grad \E_{x, y \sim p} \left[
\log p(x, y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))
\right]\\
&= \grad \E_{y \sim p} [ \E_{x \sim p(\cdot | y) } [ \log p(x, y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x)) ] ]\\
&= \grad \E_{y \sim p} [ \UU(p, y, \mathcal{S}) ].
\end{align*}

Finally, we consider $\E_{y \sim p} [ \E [ \mathbf{g} \cdot R(y) | y ] ]$ (and recall that $R(y)$ is not to be treated as a function of $\theta$):

\begin{align*}
\E_{y \sim p} [ \E [ \mathbf{g} \cdot R(y) | y ] ]
&= \E_{y \sim p} \left[ \E [ (\grad \log p(x, y)) \cdot R(y) | y ] \right]\\
&= \E_{x, y \sim p} \left[ (\grad \log p(x, y)) \cdot R(y) \right]\\
&= \grad \E_{x, y \sim p} \left[ R(y) \right]\\
&= \grad \E_{y \sim p} [ R(y) ].
\end{align*}

§.§ Proof of Theorem 3

Theorem 3. Consider an unnormalized target distribution $\tilde{\pi}(x) = Z\pi(x)$ and an inference strategy $\mathcal{S}$ targeting $\pi(x)$. Then the relative variances of the estimators $\hat{Z}(\tilde{\pi}, \mathcal{S})$ and $\check{Z}(\tilde{\pi}, \mathcal{S})$ are given by the following recursive equations:

\begin{align*}
\text{Var}_{\hat{Z}}&(\pi, \mathcal{S}) = \chi^2(\pi || \mathcal{S}.q) \, +\\
& \mathbb{E}_{x \sim \mathcal{S}.q}\left[\left(\frac{\pi(x)^2}{\mathcal{S}.q(x)^2}\right) \cdot \text{Var}_{\check{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right]\\
\text{Var}_{\check{Z}}&(\pi, \mathcal{S}) = \chi^2(\mathcal{S}.q || \pi) + \\
& \mathbb{E}_{x \sim \pi}\left[\left(\frac{\mathcal{S}.q(x)^2}{\pi(x)^2}\right) \cdot \text{Var}_{\hat{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right]
%\text{Var}\left(\frac{\hat{Z}}{Z}\right) = \sum_{i=0}^{n-1} \mathbb{E}%_{x_0, \dots, x_{i-1}}
%\frac{\hat{Z}_i^2}{Z^2} \cdot \chi^2(&p_{2\lceil \frac{i}{2}\rceil}(x_{i+1} \mid x_{i}, \dots, x_0) \,\, ||\\
%&p_{2\lfloor\frac{i}{2}\rfloor + 1}(x_{i+1} \mid x_i, \dots, x_0))
%w_i^2 \cdot d_i%
\end{align*}

When $\mathcal{S}.q$ is tractable, the second term of each sum is 0.

Proof. The proof is by induction on the level of nesting present in the inference strategy $\mathcal{S}$.

First suppose $\mathcal{S}.q$ has a tractable marginal density. Then:

* $\hat{Z}(\pi, \mathcal{S})$ is the normalized importance weight
$\frac{\pi(x)}{\mathcal{S}.q(x)}$, with $x \sim \mathcal{S}.q$. So the relative variance is:

$$\text{Var}_{\hat{Z}}({\pi}, \mathcal{S}) = \text{Var}\left(\hat{Z}({\pi}, \mathcal{S})\right) = \mathbb{E}_{x \sim \mathcal{S}.q}\left[\frac{\pi(x)^2}{\mathcal{S}.q(x)^2}\right] - \mathbb{E}_{x \sim \mathcal{S}.q}\left[\frac{\pi(x)}{\mathcal{S}.q(x)}\right]^2 = \mathbb{E}_{x \sim \mathcal{S}.q}\left[\frac{\pi(x)^2}{\mathcal{S}.q(x)^2} - 1\right] = \chi^2(\pi || \mathcal{S}.q),$$

where the third equality holds because $\pi$ is a normalized density and $\pi$ is absolutely continuous with respect to $\mathcal{S}.q$.

* $\check{Z}({\pi}, \mathcal{S})$ is the weight $\frac{\pi(x)}{\mathcal{S}.q(x)}$, with $x \sim \pi$. Then the relative variance
$$\text{Var}_{\check{Z}}({\pi}, \mathcal{S}) = \text{Var}\left(\check{Z}(\pi, \mathcal{S})^{-1}\right) = \mathbb{E}_{x \sim \pi}\left[\frac{\mathcal{S}.q(x)^2}{\pi(x)^2}\right] - \mathbb{E}_{x \sim \pi}\left[\frac{\mathcal{S}.q(x)}{\pi(x)}\right]^2 = \mathbb{E}_{x \sim \pi}\left[\frac{\mathcal{S}.q(x)^2}{\pi(x)^2} - 1\right] = \chi^2(\mathcal{S}.q || \pi),$$
where the third equality holds because $\mathcal{S}.q$ is a normalized density and is absolutely continuous with respect to $\pi$.

Now consider the inductive step. Assume that for all $x$, the theorem holds of the strategy $\mathcal{S}.\mathcal{M}(x)$ targeting $\mathcal{S}.q(\cdot \mid x)$. $\text{Var}_{\hat{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))$, for all $x$. Then:

* The $\texttt{IMPORTANCE}(\pi, \mathcal{S})$ algorithm generates $x \sim \mathcal{S}.q$. It then calls  (with $r \sim \mathcal{S}.q(\cdot \mid x)$) to obtain $w = \check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))^{-1}$, and returns $\hat{Z} = w\pi(x)$. The variance of $\hat{Z}$ is then:

\begin{align*}
\text{Var}_{\hat{Z}}(\pi, \mathcal{S}) &= \text{Var}\left(\frac{\pi(x)}{\check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))}\right)\\
&= \mathbb{E}\left[\left(\frac{\pi(x)}{\check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))}\right)^2 - 1\right]
&& (\mathbb{E}[\hat{Z}(\pi, \mathcal{S})]^2 = Z^2 = 1)\\
&= \mathbb{E}\left[\left(\frac{\pi(x)}{\mathcal{S}.q(x)} \cdot \frac{\mathcal{S}.q(x)}{\check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))}\right)^2 - 1\right]
&& \text{(divide and multiply by } \mathcal{S}.q(x))\\
&= \mathbb{E}\left[\left(\frac{\pi(x)}{\mathcal{S}.q(x)} \cdot \frac{1}{\check{Z}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))}\right)^2 - 1\right]
&&(\mathcal{S}.q(x) \text{ is the normalizing constant of } \mathcal{S}.q(\cdot, x))\\
&= \mathbb{E}\left[\left(\frac{\pi(x)}{\mathcal{S}.q(x)}\right)^2 \left(\mathbb{E}\left[{\check{Z}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))}^{-2} \bigl\vert x\right]\right) - 1\right]\\
&= \mathbb{E}\left[\left(\frac{\pi(x)}{\mathcal{S}.q(x)}\right)^2 \left(\text{Var}_{\check{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x)) + 1\right) - 1\right]
&& \text{(definition of Var}_{\check{Z}}(\cdot, \cdot))\\
&= \mathbb{E}\left[\left(\frac{\pi(x)}{\mathcal{S}.q(x)}\right)^2 \left(\text{Var}_{\check{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right) + \left(\frac{\pi(x)}{\mathcal{S}.q(x)}\right)^2 - 1\right]
&& \text{(distributing product over sum)}\\
&= \mathbb{E}\left[\left(\frac{\pi(x)}{\mathcal{S}.q(x)}\right)^2 \left(\text{Var}_{\check{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right)\right] + \chi^2(\pi || \mathcal{S}.q).\\
\end{align*}

* The argument for $\check{Z}$ is largely the same:

\begin{align*}
\text{Var}_{\check{Z}}(\pi, \mathcal{S}) &= \text{Var}\left(\frac{\hat{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))}{\pi(x)}\right)\\
&= \mathbb{E}\left[\left(\frac{\check{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))}{\pi(x)}\right)^2 - 1\right]
&& (\mathbb{E}[\check{Z}(\pi, \mathcal{S})^{-1}]^2 = Z^{-2} = 1)\\
&= \mathbb{E}\left[\left(\frac{\mathcal{S}.q(x)}{\pi(x)} \cdot \frac{\hat{Z}(\mathcal{S}.q(\cdot, x), \mathcal{S}.\mathcal{M}(x))}{\mathcal{S}.q(x)}\right)^2 - 1\right]
&& \text{(divide and multiply by } \mathcal{S}.q(x))\\
&= \mathbb{E}\left[\left(\frac{\mathcal{S}.q(x)}{\pi(x)} \cdot {\hat{Z}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))}\right)^2 - 1\right]
&&(\mathcal{S}.q(x) \text{ is the normalizing constant of } \mathcal{S}.q(\cdot, x))\\
&= \mathbb{E}\left[\left(\frac{\mathcal{S}.q(x)}{\pi(x)}\right)^2 \left(\mathbb{E}\left[{\hat{Z}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))}^{2} \bigl\vert x\right]\right) - 1\right]\\
&= \mathbb{E}\left[\left(\frac{\mathcal{S}.q(x)}{\pi(x)}\right)^2 \left(\text{Var}_{\hat{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x)) + 1\right) - 1\right]
&& \text{(definition of Var}_{\hat{Z}}(\cdot, \cdot))\\
&= \mathbb{E}\left[\left(\frac{\mathcal{S}.q(x)}{\pi(x)}\right)^2 \left(\text{Var}_{\hat{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right) + \left(\frac{\mathcal{S}.q(x)}{\pi(x)}\right)^2 - 1\right]
&& \text{(distributing product over sum)}\\
&= \mathbb{E}\left[\left(\frac{\mathcal{S}.q(x)}{\pi(x)}\right)^2 \left(\text{Var}_{\hat{Z}}(\mathcal{S}.q(\cdot \mid x), \mathcal{S}.\mathcal{M}(x))\right)\right] + \chi^2(\mathcal{S}.q || \pi).\\
\end{align*}

§.§ Proof of Theorem 4.

Theorem 4.
Consider a joint distribution $p(x, y)$ and an inference strategy $\mathcal{S}$ targeting $p(x \mid y)$. Then the following equations give the bias of ${\hat{\LL}}$ and ${\hat{\UU}}$ as estimators of $\log p(y)$:

\begin{align*}
\text{Bias}_\mathcal{L}(p, y, \mathcal{S}) =&\, -\text{KL}(\mathcal{S}.q || p(\cdot \mid y)) \\
&-\mathbb{E}_{x \sim \mathcal{S}.q}[\text{Bias}_\mathcal{U}(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
\text{Bias}_\mathcal{U}(p, y, \mathcal{S}) =&\, \text{KL}(p(\cdot \mid y) || \mathcal{S}.q)\\
&\,\,\,-\mathbb{E}_{x \sim p(\cdot \mid y)}[\text{Bias}_\mathcal{L}(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]
\end{align*}

where the second term in each equation is 0 when $\mathcal{S}.q$ has a tractable marginal density.

In the base case, where $\mathcal{S}.q$ has a tractable marginal density, the theorem states that $\log p(y) - \LL(p, y, \mathcal{S}) = KL(\mathcal{S}.q || p(\cdot \mid y))$, the familiar relationship between the standard ELBO and the KL divergence. The $\UU$ case is similar:

\begin{align*}
\text{Bias}_\UU(p, y, \mathcal{S})
&= \E_{x \sim p(\cdot \mid y)}[\log p(x, y) - \log \mathcal{S}.q(x)] - \log p(y)\\
&= \log p(y) + \E_{x \sim p(\cdot \mid y)}[\log p(x \mid y) - \log \mathcal{S}.q(x)] - \log p(y)\\
&= KL(p(\cdot \mid y) || \mathcal{S}.q).
\end{align*}

Now consider the inductive step, in which $\mathcal{S}.q$ does not have a tractable marginal density. We assume the theorem holds for $\mathcal{S}.q$ and $\mathcal{S}.\mathcal{M}(x)$. Then:

\begin{align*}
\text{Bias}_\LL(p, y, \mathcal{S})
&= \LL(p, y, \mathcal{S}) - \log p(y)\\
&= \E_{x \sim \mathcal{S}.q}[\log p(x, y) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))] - \log p(y)\\
&= \log p(y) + \E_{\sim \mathcal{S}.q}[\log p(x \mid y) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= \E_{x \sim \mathcal{S}.q}[\log p(x \mid y) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= \E_{x \sim \mathcal{S}.q}[\log p(x \mid y) - \log \mathcal{S}.q(x) + \log \mathcal{S}.q(x) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= -KL(\mathcal{S}.q || p(\cdot \mid y)) + \E_{x \sim \mathcal{S}.q}[\log \mathcal{S}.q(x) - \UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= -KL(\mathcal{S}.q || p(\cdot \mid y)) - \E_{x \sim \mathcal{S}.q}[\text{Bias}_\UU(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))].
\end{align*}

Nearly the same proof applies for $\UU$, flipping the necessary signs:

\begin{align*}
\text{Bias}_\UU(p, y, \mathcal{S})
&= \UU(p, y, \mathcal{S}) - \log p(y)\\
&= \E_{x \sim p(\cdot \mid y)}[\log p(x, y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))] - \log p(y)\\
&= \log p(y) + \E_{\sim p(\cdot \mid y)}[\log p(x \mid y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= \E_{x \sim p(\cdot \mid y)}[\log p(x \mid y) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= \E_{x \sim p(\cdot \mid y)}[\log p(x \mid y) - \log \mathcal{S}.q(x) + \log \mathcal{S}.q(x) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= KL(p(\cdot \mid y) || \mathcal{S}.q) + \E_{x \sim p(x \mid y)}[\log \mathcal{S}.q(x) - \LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))]\\
&= KL(p(\cdot \mid y) || \mathcal{S}.q) - \E_{x \sim p(\cdot \mid y)}[\text{Bias}_\LL(\mathcal{S}.q, x, \mathcal{S}.\mathcal{M}(x))].
\end{align*}

§.§ Stationarity of MCMC algorithm

In Section <ref>, we mention that RAVI can be used to run Metropolis-Hastings kernels
with proposals that have intractable densities.
Here, we present and justify the algorithm.

Let $\tilde{\pi}(x) = \int \tilde{\pi}(r, x) \text{d}r = Z \int \pi(r, x) \text{d}r$ be a possibly unnormalized target density, and let $q(x'; x) = \int q(s, x'; x) \text{d}s$ be a proposal kernel mapping previous state $x$ to new state $x'$. We note that (1) both $\tilde{\pi}$ and $q$ have intractable marginal densities, and (2) the target marginal $\tilde{\pi}(x)$ itself may be unnormalized. As is typical in pseudomarginal MCMC, even this unnormalized target density cannot be evaluated pointwise, due to the additional nuisance variables $r$.

Now suppose we have a family of inference strategies $\mathcal{S}(x)$ targeting $\pi(r \mid x)$, and a family of inference strategies $\mathcal{M}(x, x')$ targeting $q(s \mid x'; x)$. Let $x$ be a starting position for our Markov chain. We can run Algorithm [alg:alg1]1 on $\mathcal{S}$, targeting $\pi(r \mid x)$, to obtain an initial estimate $\hat{Z}_x$ of the unnormalized marginal density $\tilde{\pi}(x)$. Then Algorithm [alg:mh]5 defines a stationary MCMC kernel for the target distribution $\pi(x)$, starting at input point $x$:

Algorithm 5: RAVI Metropolis-Hastings
model $\tilde{\pi}(x) = Z \int \pi(r, x) \text{d}r$
proposal $q(x'; x) = \int q(s, x'; x) \text{d}s$
family $\mathcal{S}(x)$ of inference strategies targeting $\pi(r \mid x)$
family $\mathcal{M}(x, x')$ of inference strategies targeting $q(s \mid x'; x)$
initial position $x$ and estimate $\hat{Z}_{x}$ of $\tilde{\pi}(x)$
next position $x'$ and estimate $\hat{Z}_{x'}$ of $\tilde{\pi}(x')$
$(s, x') \sim q(s, x'; x)$
$w_{x'} \gets {\texttt{HME}}(q(\cdot, x'; x), s, \mathcal{M}(x, x'))$
$(\_, w_x) \gets \texttt{IMPORTANCE}(q(\cdot \mid x; x'), \mathcal{M}(x', x))$
$(\_, \hat{Z}_{x'}) \gets \texttt{IMPORTANCE}(\pi(\cdot \mid x'), \mathcal{S}(x'))$
$u \sim \text{Uniform}(0, 1)$
$u < \text{min}(1, \frac{\hat{Z}_{x'}}{\hat{Z}_x}w_{x'}w_x)$
$(x', \hat{Z}_{x'})$

$(x, \hat{Z}_x)$

When $q$'s marginal density is known exactly, the above algorithm recovers variants of Particle-Marginal MH [Andrieu et al., 2010], except instead of using SMC to marginalize $r$, any RAVI algorithm can be applied. When $q$'s marginal density is unavailable, however, the algorithm instead becomes a pseudo-marginal ratio algorithm [Andrieu et al., 2018], because not just $p$ but also $q$ is estimated unbiasedly. In general, it is not valid to use arbitrary unbiased estimates of $p$ and $q$, or even of $\alpha = \frac{p(x')q(x; x')}{p(x)q(x'; x)}$, within an MH algorithm. However, the added structure of the RAVI strategy ensures that the above procedure is sound.

To see why our MCMC kernel is stationary, we consider an extended target distribution. First, some notation. For an inference strategy $\mathcal{S}$ targeting $\pi(x)$, write $v_\mathcal{S}$ for the complete set of auxiliary variables in the strategy: if $\mathcal{S}.q$ has a tractable marginal density, then $v_\mathcal{S} = \emptyset$, and otherwise, if $\mathcal{S}.q(x) = \int \mathcal{S}.q(r, x) \text{d}r$, then $v_\mathcal{S}$ is defined recursively as $\{r\} \cup v_{\mathcal{S}.\mathcal{M}}$. Calling $\texttt{IMPORTANCE}$ on $\mathcal{S}$ yields a joint distribution over these auxiliary variables and $x$, which we denote as $p^{\mathcal{S}}_{\texttt{IMP}}(v_\mathcal{S}, x)$. Calling $\texttt{HME}$ on $\mathcal{S}$ and a particular sample $x$ yields a distribution over just $v_\mathcal{S}$, which we denote $p^\mathcal{S}_{\texttt{HME}}(v_\mathcal{S}; x)$.
When $x \sim \pi$ and $v_\mathcal{S} \sim p^\mathcal{S}_\texttt{HME}(v_\mathcal{S}; x)$, the ratio $\frac{p^\mathcal{S}_\texttt{IMP}(v_\mathcal{S}, x)}{\tilde{\pi}(x)p^\mathcal{S}_\texttt{HME}(v_\mathcal{S}; x)}$ is the weight $\check{Z}(\tilde{\pi}, \mathcal{S})^{-1}$ returned by , and similarly, when $(v_\mathcal{S}, x) \sim p^\mathcal{S}_\texttt{IMP}$, the ratio $\frac{\tilde{\pi}(x)p^\mathcal{S}_\texttt{HME}(v_\mathcal{S}; x)}{p^\mathcal{S}_\texttt{IMP}(v_\mathcal{S}, x)}$ is the weight $\hat{Z}(\tilde{\pi}, \mathcal{S})$ returned by $\texttt{IMPORTANCE}$.

Using this notation, we can extend the target distribution $\tilde{\pi}(x)$ to one over $(x, s, x', s', v_{\mathcal{S}(x)}, v_{\mathcal{M}(x, x')}, v_{\mathcal{M}(x', x)})$ that admits $\tilde{\pi}(x)$ as a marginal:

$$\tilde{\pi}(r, x, s, x', s', v_{\mathcal{S}(x)}, v_{\mathcal{M}(x, x')}, v_{\mathcal{M}(x', x)}) = \tilde{\pi}(r, x) \cdot p^{\mathcal{S}(x)}_\texttt{HME}(v_{\mathcal{S}(x)}; r) \cdot q(s, x'; x) \cdot p^{\mathcal{M}(x, x')}_\texttt{HME}(v_{\mathcal{M}(x, x')}; s) \cdot p^{\mathcal{M}(x', x)}_\texttt{IMP}(v_{\mathcal{M}(x', x)}, s')$$

Our algorithm can be understood as sequencing two stationary kernels for this extended target. The first (implemented by lines 1-3) is a blocked Gibbs update on the variables $(s, x', s', v_{\mathcal{M}(x, x')}, v_{\mathcal{M}(x', x)})$, conditioned on everything else. Lines 1-3 sample exactly from the conditional distribution of these variables. The second is a Metropolis-Hastings proposal that simultaneously: (i) swaps $x$ with $x'$ (the `main' proposed update), (ii) swaps $(s, v_{\mathcal{M}(x, x')})$ with $(s', v_{\mathcal{M}(x', x)})$, and (iii) proposes an update to $r$ and to $v_{\mathcal{S}(x)}$ from $p_\texttt{IMP}^{\mathcal{S}(x')}$. The usual Metropolis-Hastings acceptance probability for this kernel, computed on the extended state space, is precisely the formula in Line 6.

One consequence of this justification is that the same family $\mathcal{S}$ of inference strategies for $\pi$ must be used at each iteration. The family $\mathcal{M}$ can be freely switched out (as can $q$), however, to develop a cycle of kernels that use different proposal distributions.

§ FURTHER EXAMPLES

This appendix lists examples of popular Monte Carlo and variational inference algorithms, and explains how they can be viewed as inference strategies. In addition, some of these algorithms can be viewed as inference strategy combinators, because they feature user-chosen proposal distributions or variational families that can themselves be instantiated with inference strategies.[This `combinator' viewpoint evokes earlier work by [Ś]cibior et al., 2017] and [Stites et al., 2021]. For example, [Stites et al., 2021] introduce combinators for creating properly weighted samplers compositionally, with parameters that can be optimized using standard or nested variational objectives. Some of their combinators have equivalents in this section, e.g. their  combinator is similar to the construction we present for Nested Importance Sampling in Section <ref>. However: (1) the fundamental compositional operation in RAVI, of combining a posterior approximation with a meta-posterior approximation, cannot be achieved using their combinators; (2) as such, some of the algorithms that RAVI covers cannot be constructed using their combinators; and (3) their combinators produce properly weighted samplers, which contain `less information' than inference strategies: an inference strategy can be used, e.g., as a proposal distribution in Metropolis-Hastings, whereas properly weighted samplers cannot in general be used this way.

§.§ $N$-particle Importance Sampling

RAVI Inference Strategy: $N$-particle Importance Sampling
sirsir($\tilde{\pi}, q, N$).qsirmsir($\tilde{\pi}, q, N$).M($x$).q

infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

latent variable $x$
particles $x_{1:N}$, chosen particle index $j$
$i \in 1, \dots, N$
$x_i \sim q$
$w_i \gets \frac{\tilde{\pi}(x_i)}{q(x_i)}$

$j \sim \text{Discrete}(w_{1:N})$

metaalgMeta-Posterior Approx.

particles $x_{1:N}$, chosen particle index $j$
$j \sim \text{Uniform}(1, N)$
$x_j \gets x$
$i \in 1,\dots,j-1,j+1,\dots,N$
$x_i \sim q$

$(x_{1:N}, j)$

Consider the $N$-particle importance sampling estimator
$$\hat{Z} = \frac{1}{N} \sum_{i=1}^N \frac{\tilde{\pi}(x_i)}{q(x_i)},\text{ for }x_i \sim q.$$
The same estimator can be recovered as a one-particle  estimate, by applying Alg. [alg:alg1]1 to the [infstrat:sir] inference strategy.

The proposal $\mathcal{S}.q$ generates $N$ particles $x_{1:N}$, and selects an index $j$ from a discrete distribution on $1, \dots, N$, with weights proportional to $w_i = \tilde{\pi}(x_i) / q(x_i)$. The meta-proposal is responsible for inferring $j$ and the complete set of particles $x_{1:M}$, given the chosen particle $x$. It uses the conditional SIR algorithm [Andrieu et al., 2010] to do so, proposing $j$ uniformly in $\{1, \dots, N\}$, and generating values for the un-chosen particles $x_{-j}$ from $q$.

RAVI Inference Strategy: $N$-particle IS with RAVI strategy $\mathcal{S}$
ravisirravi-sir($\tilde{\pi}, \mathcal{S}, N$).qravisirmravi-sir($\tilde{\pi}, \mathcal{S}, N$).M($x$).q

infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

latent variable $x$
particles $x_{1:N}$, aux. proposal variables $v_\mathcal{S}^{1:N}$, chosen particle index $j$
$i \in 1, \dots, N$
$x_i, w_i \sim \texttt{IMPORTANCE}(\tilde{\pi}, \mathcal{S})$ w. aux. vars $v_\mathcal{S}^i$

$j \sim \text{Discrete}(w_{1:N})$

metaalgMeta-Posterior Approx.

particles $x_{1:N}$, aux. proposal variables $v_\mathcal{S}^{1:N}$, chosen particle index $j$
$j \sim \text{Uniform}(1, N)$
$x_j \gets x$
$\_ \sim \texttt{HME}(\tilde{\pi}, x_j, \mathcal{S})$ w. aux. vars $v_\mathcal{S}^j$
$i \in 1,\dots,j-1,j+1,\dots,N$
$\_, x_i \sim q$ w. aux. vars $v_{\mathcal{S}}^i$

$(v_\mathcal{S}^{1:N}, x_{1:N}, j)$

This is a suboptimal choice of $\mathcal{S}.\mathcal{M}(x).q$; lower-variance estimates $\hat{Z}$ can be obtained by improving meta-inference, either by incorporating problem-specific domain knowledge or via learning. However, in many cases, improved meta-inference may not be worth the computation required; it remains to be seen whether techniques such as amortized learning can be applied to deliver accuracy gains at low computational cost.

Instantiating the proposal $q$ as its own inference strategy.
The above assumes that $q$ has a tractable marginal density. When it doesn't, the inner importance sampling loop can use a RAVI inference strategy $\mathcal{S}$ instead of a tractable proposal $q$. This modification is presented in the higher-order inference strategy [infstrat:ravisir].
One way to think about this construction is as a way to improve any existing inference strategy $\mathcal{S}$ by `adding replicates.' The resulting estimator of $Z$ is the mean of $N$ independent $\hat{Z}$ estimates from the original inference strategy.

§.§ Importance-Weighted Autoencoders

The importance-weighted auto-encoder arises by considering the same inference strategy as in Section <ref>, but as a variational inference procedure (Alg. 3) rather than a Monte Carlo procedure.

Because $\texttt{sir}(\tilde{\pi}, q, N).q$ of this inference strategy corresponds to $N$-particle sampling importance-resampling (SIR), it has been argued that IWAE is in fact `vanilla' variational inference, but with a variational family that uses SIR to more closely approximate the posterior [Bachman and Precup, 2015]. However, [Cremer et al., 2017] show that deriving the ELBO for that variational family gives rise to a different objective, and that IWAE gives a looser lower bound on $\log Z$ than this idealized (but generally intractable) objective.

In the RAVI framework, these two objectives arise from different inference strategies, which share the same $\mathcal{S}.q$ (SIR in both cases), but use different meta-inference $\mathcal{S}.\mathcal{M}$. IWAE uses the simple conditional SIR meta-inference introduced in Section <ref>, whereas [Cremer et al., 2017]'s idealized objective can be derived by using the optimal choice of $\mathcal{S}.\mathcal{M}(x).q(j, x_{1:N})$—the exact posterior of the SIR procedure. The looser bound obtained by IWAE can be seen as a result of its $\mathcal{S}.\mathcal{M}$ performing poorer meta-inference: inference about the auxiliary variables of the SIR inference algorithm used in $\mathcal{S}.q$.

§.§ $N$-particle Sequential Monte Carlo

RAVI Inference Strategy: $N$-particle SMC w. RAVI strategies
ravismcsmc($\tilde{\pi}_{1:T}, \mathcal{S}, K_{2:T}, L_{2:T}, N$).qravismcmsmc($\tilde{\pi}_{1:T}, \mathcal{S}, K_{2:T}, L_{2:T}, N$).M($x$).q

infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

latent variable $x$ targeting $\tilde{\pi}_T$
particles $x^{1:T}_{1:N}$, aux. proposal variables $v_\mathcal{S}^{1:N}$, aux. $K$ vars $v_{K_{2:T}}^{1:N}$, aux. $L$ vars $v_{L_{2:T}}^{1:N}$, ancestor variables $a^{1:T-1}_{1:N}$, final chosen particle index $j$
$i \in 1, \dots, N$
$x^1_i, w^1_i \sim \texttt{IMPORTANCE}(\tilde{\pi}_1, \mathcal{S})$ w. aux. vars $v_\mathcal{S}^i$

$t \in 2, \dots, T$
$i \in 1, \dots, N$
$a_i^{t-1} \sim \text{Discrete}(w^{t-1}_{1:N})$
$x_i^t, \hat{w} \sim \texttt{IMPORTANCE}(\tilde{\pi}_t, K_t(x^{t-1}_{a^{t-1}_i}))$ w. aux. vars $v_{K_t}^i$
$\check{w} \sim \texttt{HME}(\tilde{\pi}_{t-1}, x^{t-1}_{a^{t-1}_i}, L_t(x_i^t))$ w. aux. vars $v_{L_t}^i$
$w_i^{t} \gets \hat{w} \cdot \check{w}$

$j \sim \text{Discrete}(w^T_{1:N})$

metaalgMeta-Posterior Approx.

particles $x^{1:T}_{1:N}$, aux. proposal variables $v_\mathcal{S}^{1:N}$, aux. $K$ vars $v_{K_{2:T}}^{1:N}$, aux. $L$ vars $v_{L_{2:T}}^{1:N}$, ancestor variables $a^{1:T-1}_{1:N}$, final chosen particle index $j$
$j \sim \text{Uniform}(1, N)$
$x_j^T, b_T \gets x, j$
$t \in T, \dots, 2$
$a_{b_t}^{t-1} \sim \text{Uniform}(1, N)$
$b_{t-1} \gets a_{b_t}^{t-1}$
$x_{b_{t-1}}^{t-1}, \check{w} \sim \texttt{IMPORTANCE}(\tilde{\pi}_{t-1}, L_t(x_{b_t}^t))$ w. aux. vars $v_{L_t}^{b_t}$
$\hat{w} \sim \texttt{HME}(\tilde{\pi}_{t}, x_{b_t}^t, K_t(x_{b_{t-1}}^{t-1}))$ w. aux. vars $v_{K_t}^{b_t}$
$w_{b_t}^t \gets (\hat{w} \cdot \check{w})^{-1}$

$w^1_{b_1} \sim \texttt{HME}(\tilde{\pi}_1, x_{b_1}^1, \mathcal{S})$ w. aux. vars $v_\mathcal{S}^{b_1}$
$i \in 1, \dots, b_1-1, b_1+1, \dots, N$
$x_{i}^1, w_i^1 \sim \texttt{IMPORTANCE}(\tilde{\pi}_1, \mathcal{S})$ w. aux. vars $v_\mathcal{S}^i$

$t \in 2, \dots, T$
$i \in 1, \dots, b_t-1, b_t+1, \dots, N$
$a_i^{t-1} \sim \text{Discrete}(w^{t-1}_{1:N})$
$x_i^t, \hat{w} \sim \texttt{IMPORTANCE}(\tilde{\pi}_t, K_t(x^{t-1}_{a^{t-1}_i}))$ w. aux. vars $v_{K_t}^i$
$\check{w} \sim \texttt{HME}(\tilde{\pi}_{t-1}, x^{t-1}_{a^{t-1}_i}, L_t(x_i^t))$ w. aux. vars $v_{L_t}^i$
$w_i^{t} \gets \hat{w} \cdot \check{w}$

$(x_{1:N}^{1:T}, v_\mathcal{S}^{1:N}, v_{K_{2:T}}^{1:N}, v_{L_{2:T}}^{1:N}, a_{1:N}^{1:T-1},j)$

The sequential Monte Carlo family of algorithms [Chopin and Papaspiliopoulos, 2020, Del Moral et al., 2006] evolve a population of weighted particles to approximate a sequence of target distributions. SMC can be viewed as standard importance sampling, with an inference strategy in which $\mathcal{S}.q$ is the sampling distribution for SMC, and $\mathcal{S}.\mathcal{M}(x)$ is the conditional SMC algorithm [Andrieu et al., 2010].

Standard SMC is parameterized by:

* A sequence $\tilde{\pi}_{1:T}$ of intermediate target distributions, with $\tilde{\pi}_T = \tilde{\pi}$ the ultimate target;

* An initial proposal $q(x_1)$;

* A sequence $K_t(x_{t-1} \rightarrow x_t)$ of proposal kernels for $t=2, \dots, T$; and

* A sequence $L_t(x_{t} \rightarrow x_{t-1})$ of backward kernels for $t=2, \dots, T$.

Here, we show a version of SMC (the inference strategy [infstrat:smc]) that behaves as a `higher-order inference strategy,' or `inference strategy combinator': it allows for an initial proposal, proposal kernels, and backward kernels that do not have tractable marginal densities. Our version is parameterized by:

* A sequence $\tilde{\pi}_{1:T}$ of intermediate target distributions, with $\tilde{\pi}_T = \tilde{\pi}$ the ultimate target;

* An initial proposal $\mathcal{S}$ (a RAVI strategy);

* A sequence of inference strategy families $K_t(x_{t-1})$ parameterized by $x_{t-1}$, for $t=2, \dots, T$, targeting $\tilde{\pi}_t$; and

* A sequence of inference strategy families $L_t(x_{t})$ of backward kernels, parameterized by $x_t$, for $t=2, \dots, T$.

The posterior approximation $\mathcal{S}.q$ runs a version of SMC that uses $\texttt{HME}$ and $\texttt{IMPORTANCE}$ to compute weights. The meta-posterior approximation $\mathcal{S}.\mathcal{M}(x).q$ runs a similarly modified version of conditional SMC [Andrieu et al., 2010]. When  is run on the  inference strategy, the final weight $\hat{Z}$ is the SMC marginal likelihood esitmate, the product of the averages of the weights from each time step.

It is possible to adapt this strategy to use adaptive resampling and rejuvenation. (Rejuvenation moves do not actually require modification: can be incorporated by including them as explicit $(K, L)$ pairs, where $L$ is the time-reversal of an MCMC kernel $K$.) However, we are not aware of a way to justify the adaptive choice of rejuvenation kernel.

§.§ Variational Sequential Monte Carlo

The Variational Sequential Monte Carlo [Naesseth et al., 2018] objective corresponds to Alg. 3, with the same RAVI inference strategy as in Appendix <ref>. However, the default gradient estimator from Alg. 3 will have high variance. Naesseth et al., 2018 recommend using a biased estimator of the gradient, that uses reparameterization where possible and discards the score function terms arising from resampling steps.

§.§ Annealed Importance Sampling

RAVI Inference Strategy: Annealed Importance Sampling
aisais($\tilde{\pi}_{1:T}, \mathcal{S}, K_{2:T}$).qaismais($\tilde{\pi}_{1:T}, \mathcal{S}, K_{2:T}$).M($x$).q

infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

latent variable $x$ targeting $\tilde{\pi}_T$
$x^{1:T}$, aux. vars $v_\mathcal{S}$ of initial proposal
$x_1, \_ \sim \texttt{IMPORTANCE}(\tilde{\pi}_1, \mathcal{S})$ w. aux. vars $v_\mathcal{S}$
$t \in 2, \dots, T$
$x_t \sim K_t(x_{t-1} \rightarrow \cdot)$

metaalgMeta-Posterior Approx.

$x^{1:T}$, aux. vars $v_\mathcal{S}$ of initial proposal
$x_T \gets x$
$t \in T, \dots, 2$
$x_{t-1} \sim \tilde{K}_t(x_t \rightarrow \cdot)$*$\tilde{K}_t$ is time reversal of $K_t$

$\_ \sim \texttt{HME}(\tilde{\pi}_1, x_1, \mathcal{S})$ w. aux. vars $v_\mathcal{S}$

$(x_{1:T}, v_\mathcal{S})$

In annealed importance sampling, the practitioner chooses a sequence of unnormalized target distributions $\tilde{\pi}_{1:T}$, where $\pi_T$ is the posterior distribution of interest. Typically $\pi_1$ is chosen to be a distribution that is easy to approximate with a proposal $q$, and each $\pi_i$ is slightly closer to the true target $\pi_T$ than the last. The user also chooses a sequence of kernels $K_t(x_{t-1} \rightarrow x_t)$, where $K_t$ is stationary for $\pi_{t-1}$. The algorithm begins by sampling an initial point $x_1 \sim q$, transforming it through the sequence of kernels to obtain $x_2, \dots, x_T$, and returning $x_T$ as the inferred value of $x$. The associated weight is
$$\hat{Z} = \frac{\tilde{\pi}_1(x_1) \cdot \dots \cdot \tilde{\pi}_T(x_T)}{q(x_1) \cdot \tilde{\pi}_1(x_2) \cdot \dots \cdot \tilde{\pi}_{T-1}(x_T)}.$$

This procedure corresponds to running Alg. [alg:alg1]1 on the [infstrat:ais] inference strategy. The inference process runs the kernels $K_t$ forward, whereas the meta-inference process runs their time reversals backward: $\tilde{K}_t(x_t \rightarrow x_{t-1}) \propto \pi_t(x_{t-1}) \cdot K_t(x_{t-1} \rightarrow x_t)$.

Note that if $K$ is a stationary kernel for $\pi_i$, so is $K^m$ for any natural number $m$. With sufficient computation (increasing $m$), we can ensure that the AIS top-level proposal $\texttt{ais}(\dots).q$ is arbitrarily close to the target posterior $\pi_T$. However, doing so will not necessarily lead to lower-variance weights: RAVI makes clear that it is also necessary to consider the quality of meta-inference.

Consider the job of $\tilde{K}_T$, which in the context of the meta-posterior approximation $\texttt{ais}.\mathcal{M}(x)$ is supposed to infer $x_{T-1}$ from $x_{T}$. $\tilde{K}_T$  is the exact meta-posterior of $x_{T-1}$ given $x_T$ assuming that, in the forward direction, $x_{T-1}$ was distributed according to $\pi_{T-1}$. However, in the forward direction, if each $K_t$ is run sufficiently many times to ensure mixing at each step, $x_{T-1}$ will in fact be distributed according to $\pi_{T-2}$. This gap—between the optimal meta-inference kernels and the actual $\tilde{K}$ kernels—is partly responsible for the variance of the AIS estimator, and can be mitigated by using a finer annealing schedule that brings successive target distributions closer together.  It could also be mitigated by learning a better reverse annealing chain.

§.§ Nested Sequential Monte Carlo

We first consider Nested Importance Sampling. As in RAVI, Nested Importance Sampling is concerned with importance sampling when the proposal distribution $q$  cannot be tractably evaluated. But RAVI and NIS take different approaches:

* RAVI assumes $q$ can be simulated, but that the (normalized) density cannot be evaluated. RAVI generates proposals exactly distributed according to the user's desired proposal $\mathcal{S}.q$, and generates approximations to the ideal importance weights.

* NIS does not assume $q$ can be simulated, but does assume that its unnormalized density $\tilde{q}$ is available. As such, proposals are not simulated from $q$, but rather from a Sampling/Importance-Resampling (SIR) approximation to $q$.

The NIS procedure with an intractable proposal $q$ corresponds exactly to a special case of the RAVI algorithm, with the RAVI proposal $\mathcal{S}.q$ set not to $q$ but rather to an SIR sampling distribution targeting $q$ using some tractable proposal $h$. Compare:

* Ordinary SIR targeting $\tilde{\pi}$ with proposal $h$: recovered by running $\texttt{IMPORTANCE}(\tilde{\pi}, \texttt{sir}(\tilde{\pi}, h, N))$ (see Section <ref> for  inference strategy).

* Nested IS targeting $\tilde{\pi}$ with unnormalized proposal density $\tilde{q}$, approximated using SIR with $h$ as a proposal: recovered by running $\texttt{IMPORTANCE}(\tilde{\pi}, \texttt{sir}(\tilde{q}, h, N))$.

That is, under the RAVI perspective, the only difference between ordinary SIR using $h$ and nested IS is that the ideal proposal density $\tilde{q}$ (rather than the target density $\tilde{\pi}$) is used to make the resampling decision about the particles generated by $h$ (the index $j$ in the listing for ).

More generally, Naesseth et al., 2015 consider procedures other than SIR for approximating ${q}$, arguing that any properly weighted sampler for the intractable proposal $q$ will do. If we let $\mathcal{H}$ be a RAVI inference strategy representing the properly weighted sampler for the intractable proposal $q$ (with unnormalized density $\tilde{q}$), then the Nested IS procedure that uses this properly weighted proposal to perform inference in $\tilde{\pi}$ is $\texttt{IMPORTANCE}(\tilde{\pi}, \texttt{ravi-sir}(\tilde{q}, \mathcal{H}, 1))$ (see [infstrat:ravisir] in Section <ref>).

Nested SMC is similar, performing Nested IS at each iteration of SMC. To recover this algorithm using RAVI, we use the [infstrat:smc] inference strategy, but for the proposals $K_t(x_{t-1})$ (which, as described in Section <ref>, can be instantiated with inference strategies), we use [infstrat:ravisir] targeting the desired but intractable proposal.

§.§ SMC$^2$

RAVI Inference Strategy: SMC$^2$
smcsqsmc$^2$($p, q_1, q, M, N$).qsmcsqmsmc$^2$($p, q_1, q, M, N$).M($\theta, x_{1:T}$).q

infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

parameters $\theta$, sequence $x_{1:T}$
inner SMC vars $v_\texttt{smc}^T$ of chosen SMC$^2$ particle, other SMC$^2$ vars $v$
the targets $\tilde{\pi}_t$ depend on $M$, $p$, $q_1$, and $q$
$(\theta, x_{1:T}, v_\texttt{smc}^T), \_ \sim \texttt{IMPORTANCE}(\tilde{\pi}_T, \texttt{smc}(\tilde{\pi}_{1:T}, K_{2:T}^2, L_{2:T}^2, N))$ w. aux. vars $v$
$\theta, x_{1:T}$

metaalgMeta-Posterior Approx.

inner SMC vars $v_\texttt{smc}^T$ of chosen SMC$^2$ particle, other SMC$^2$ vars $v$
$\_ \sim \texttt{HME}(p_T^\theta, x_{1:T}, \texttt{smc}(p^\theta_{1:T}, q_1, K_{2:T}, L_{2:T}, M))$ w. aux. vars $v_\texttt{smc}^T$
$\_ \sim \texttt{HME}(\tilde{\pi}_T, (\theta, x_{1:T}, v_\texttt{smc}^T), \texttt{smc}(\tilde{\pi}_{1:T}, K_{2:T}^2, L_{2:T}^2, N))$ w. aux. vars $v$

$(v_\texttt{smc}^T, v)$

Suppose we are working with a state-space model $p(\theta) \prod_{i=1}^T p(x_i \mid x_{1:i}, \theta) p(y_i \mid x_i, \theta)$. For a fixed $\theta$, an SMC algorithm could be used to target the successive posteriors $p_t^\theta(x_{1:t}) = p(x_{1:t} \mid y_{1:t}, \theta)$, with proposal kernels $K_t(x^{t-1}_{1:t-1} \rightarrow x^t_{1:t}) = \delta_{x^{t-1}_{1:t-1}}(x^t_{1:t-1})q(x^t_t; x^t_{1:t-1}, y_{1:t}, \theta)$ (for some choice of $q$) and deterministic backward kernels $L_t(x^t_{1:t} \rightarrow x^{t-1}_{1:t-1}) = \delta_{x^t_{1:t-1}}(x^{t-1}_{1:t-1})$. The RAVI strategy implementing that SMC algorithm is $\texttt{smc}(p_{1:T}, q_1, K_{2:T}, L_{2:T}, N)$, where $q_1(x_1; \theta)$ is a proposal for an initial $x_1$ and $N$ is the number of particles.

If we also wish to infer $\theta$, we can instead use the SMC$^2$ algorithm [Chopin et al., 2013]. We define extended targets

$$\pi_t(\theta, x_{1:t}, v_\texttt{smc}^t) = p(\theta \mid y_{1:t}) p(x_{1:t} \mid y_{1:t}, \theta) p_{\texttt{HME}}^{\texttt{smc}(p_{1:t}^\theta, q_1, K_{2:t}, L_{2:t}, N)}(v_\texttt{smc}^t; x_{1:t}),$$

which are defined over not only $\theta$ and $x_{1:t}$ but also all the auxiliary variables $v_\texttt{smc}^t$ used during steps 1 through $t$ of SMC. The variables $v_\texttt{smc}$ and the $p_\texttt{HME}$ distribution over them are as defined in Appendix <ref>. We write $\tilde{\pi}_t$ for the unnormalized versions of these targets, with normalizing constant $p(y_{1:t})$.

The SMC$^2$ algorithm targets this sequence of extended posteriors. We write $K_t^2$ for the forward kernels used  by this outer SMC algorithm. The kernel $K_t^2$ extends the SMC state variables $v_\texttt{SMC}^{t-1}$ to new state variables $v_\texttt{SMC}^{t}$ by running the particle filter forward one step, resampling the chosen trajectory index $j$ based on the new weights for time step $t$, and updating $x_{1:t}$ to match the $j^\text{th}$ trajectory. The corresponding backward kernel $L_t^2$ deletes the $t^\text{th}$ step of the particle deterministically, then reproposes $j$ based on the step $t-1$ weights, setting $x_{1:t-1}$ to match the $j^\text{th}$ trajectory.

The SMC$^2$ algorithm corresponds to the RAVI strategy [infstrat:smcsq]. Running the other SMC yields an approximate sample from $\tilde{\pi}_T$, which includes auxiliary variables $v_\texttt{smc}^T$. Meta-inference runs two rounds of conditional SMC: first, to recover the inner layer of SMC's variables $v_\texttt{smc}^T$ for the chosen outer-layer particle, and second, to recover the outer layer of SMC's auxiliary variables $v$.
As discussed by Chopin et al., 2013, particle MCMC rejuvenation moves can also be included; to justify using RAVI, we would insert these kernels as additional proposals within the sequence $K_{2:T}^2$.

§.§ Amortized Rejection Sampling

Consider a generative model $p(K, x_{1:K+1}, y)$ where the latent variables $x_{1:K+1}$ to be marginalized or inferred represent the trace of a rejection sampling loop, with sampling distribution $h(x)$ and predicate $\mathcal{A}(x)$ determining acceptance:

$$p(K, x_{1:K+1}, y) = \prod_{i=1}^{K} \left[h(x_i)(1-\mathcal{A}(x_i))\right] h(x_{K+1})\mathcal{A}(x_{K+1}) p(y \mid x_{K+1})$$

RAVI Inference Strategy: Amortized Rejection Sampling
amrejamrej($h, q, \mathcal{A}, N, M$).qamrejmamrej($h, q, \mathcal{A}, N, M$).M($K, x_{1:K}$).q
amrejmmamrej($h, q, \mathcal{A}, N, M$).M($K, x_{1:K}$).M($K', x'_{1:K'}, (K''_i, x''^i_{1:K''_i})_{i=1:M}, j$).q
amrejmmmamrej($h, q, \mathcal{A}, N, M$).M($K, x_{1:K}$).M($K', x'_{1:K'}, (K''_i, x''^i_{1:K''_i})_{i=1:M}, j$).M($z_{K'+1}$).q
infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

number $K$ of rejected samples, rejected samples $x_{1:K}$, accepted sample $x_{K+1}$
rejection loops $(K', x'_{1:K'})$ and $(K''_i, x''^i_{1:K''_i})_{i=1:M}$, index $j$
$K' \gets 0$
$x'_{1} \sim q$
$\mathcal{A}(x'_{K'+1}) \neq 1$
$K' \gets K' + 1$
$x'_{K'+1} \sim q$

$x_{K+1} \gets x'_{K'+1}$
$i \in 1, \dots, M$
$K''_i \gets 0$
$x''^i_{1} \sim h$
$\mathcal{A}(x''^i_{K''_i+1}) \neq 1$
$K''_i \gets K''_i + 1$
$x''^i_{K''_i+1} \sim h$

$j \sim \text{Discrete}(K''_{1:M})$
$K \sim \text{Uniform}(0, K''_j)$
$x_{1:K} \gets x''^j_{1:K}$
$(K, x_{1:K}, x_{K+1})$

metaalgMeta-Posterior Approx.

rejection loops $(K', x'_{1:K'})$ and $(K''_i, x''^i_{1:K''_i})_{i=1:M}$, index $j$
superfluous accepted sample $z_{K'+1}$
$j \sim \text{Uniform}(1, M)$
$i \in 1, \dots, j-1, j+1, \dots, M$
$K''_i \gets 0$
$x''^i_{1} \sim h$
$\mathcal{A}(x''^i_{K''_i+1}) \neq 1$
$K''_i \gets K''_i + 1$
$x''^i_{K''_i+1} \sim h$

$K''_j \gets K$
$x''^j_{K+1} \sim h$
$\mathcal{A}(x''^j_{K+K''_j+1}) \neq 1$
$K''_j \gets K''_j + 1$
$x''^j_{K+K''_j+1} \sim h$

$K' \gets 0$
$z_{1} \sim q$
$\mathcal{A}(z_{K'+1}) \neq 1$
$x'_{K'+1} \gets z_{K'+1}$
$K' \gets K' + 1$
$z_{K'+1} \sim q$

$(K', x'_{1:K'}, (K''_i, x''^i_{1:K''_i})_{i=1:M}, j)$

metaalgMeta-Meta-Posterior Approx.

superfluous accepted sample $z_{K'+1}$
index $l$, unchosen particles $z_{-l}$
$i \in 1, \dots, N$
$z_i \sim q$

$l \sim \text{Uniform}(\{i \mid \mathcal{A}(z_i)\})$

metaalgMeta-Meta-Meta-Posterior Approx.

index $l$, unchosen particles $z_{-l}$
$l \sim \text{Uniform}(1, N)$
$i \in 1, \dots, l-1,l+1,\dots, N$
$z_i \sim q$

$(z_1, \dots, z_{l-1}, z_{l+1}, \dots, z_N)$

Here, the $x_i$ are drawn independently from a distribution $h$, until some predicate $\mathcal{A}$ holds of the most recent particle, at which point the loop stops. The observation $y$ depends on the final sample $x_{K+1}$, but not the earlier, rejected samples $x_{1:K}$ or the number of rejected samples $K$.
Naderiparizi et al., 2019 proposed a technique called Amortized Rejection Sampling for performing inference in this model. The technique corresponds to the rather involved RAVI strategy , which has parameters $N$ and $M$ that can be used to trade accuracy for computational cost.

The idea behind the top-level, intractable posterior approximation $\texttt{amrej}(h, q, \mathcal{A}, N, M).q$ is to:

* use the observation $y$ to intelligently guess the accepted particle $x_{K+1}$, using a learned proposal $q$. (For example, $q$ may be parameterized by a neural network that accepts $y$ as input.) To satisfy the constraint that $x_{K+1}$ satisfies $\mathcal{A}$, however, it is necessary to run $q$ within a rejection sampling loop, generating auxiliary variables $x'_{1:K'}$, where $K'$ is the number of rejected $q$-samples. (We could try directly using $x'_{1:K'}$ as our proposal for $x_{1:K}$, the rejected samples from the model. But $q$'s goal is to propose $x_{K+1}$ in a data-driven way, influenced by the observation $y$, and the rejected samples $x_{1:K}$ from the model have no connection to the data—so, using samples from $q$ as proposals for the rejected model samples would result in a poor approximation.)

* use rejection sampling from the prior $h$ to infer the rejected samples $x_{1:K}$. We run $M$ independent rejection sampling loops, randomly choose one with probability proportional to its length, and then randomly choose a prefix of the chosen loop as our proposal for $x_{1:K}$.

The meta-posterior approximation must solve two new challenges: recovering the rejected $q$ samples $x'_{1:K'}$ from the posterior approximation, and recovering the many unused rejection loops (and the suffix of the chosen rejection loop) from the second step of the posterior approximation (the $x''$ variables). The latter of these tasks is simple enough: we can generate $M-1$ rejection loops from scratch for the un-chosen loops, and a further rejection loop from scratch to use as the suffix of the chosen loop. The first task is more complex: we run a new rejection loop using $q$ as a proposal, and discard the final accepted sample. Meta-meta-inference must infer this discarded accepted sample, for which it uses SIR with $N$ particles. The final layer, the Meta-Meta-Meta-Posterior Approximation, uses conditional SIR.

The meta-meta-posterior is not absolutely continuous with respect to its approximation (it is possible that the approximation generates $N$ $z$-values that all fail to satisfy the predicate, in which case $z_l$ is not in the support of the meta-meta-posterior). As such, this is an example of a wide inference strategy (Appendix <ref>).

§.§ Hamiltonian Variational Inference

RAVI Inference Strategy: Hamiltonian Variational Inference
hamhamvi($q_0, q_v, r_v, \text{LF}$).qhammhamvi($q_0, q_v, r_v, \text{LF}$).M($x$).q
hammmhamvi($q_0, q_v, r_v, \text{LF}$).M($x$).M($x_0, v$).q
infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

latent variable $x$
initial position $x_0$, momentum $v$
$x_0 \sim q_0$
$v \sim q_v$
$(x, v') \gets \text{LF}(x_0, v)$

metaalgMeta-Posterior Approx.

initial position $x_0$, momentum $v$
negated final momentum $v'_{-}$
$v'_{-} \sim r_v(\cdot; x)$
$(x_0, v_{-}) \gets \text{LF}(x, v'_{-})$
$(x_0, -v_{-})$

metaalgMeta-Posterior Approx.

negated final momentum $v'_{-}$
$(\_, v') \gets \text{LF}(x_0, v)$

Hamiltonian Variational Inference [Salimans et al., 2015] is a hybrid of Hamiltonian Monte Carlo and variational inference. It is a special case of Markov Chain Variational Inference (see Section <ref> and Section <ref> for detailed discussion, and [infstrat:mcvi] for the RAVI implementation). The algorithm specializes the Markov Chain Variational Inference procedure for use with a Hamiltonian Monte Carlo kernel.

We present the specialized strategy as [infstrat:ham]. It accepts as input:

* a distribution $q_0$ from which to propose an initial point;

* a momentum distribution $q_v$ from which momenta $v$ are proposed at each iteration;

* a proposal distribution $r_v(\cdot; x)$ over momenta; and

* a leapfrog integrator $\texttt{LF}$ that runs Hamiltonian dynamics on an initial position and momentum (we think of both the number of leapfrog steps $L$ and the Hamiltonian $H$ being targeted as part of the  object provided to ).

Given these inputs, the top-level posterior approximation runs an iteration of HMC from a randomly initialized location $x_0$. The meta-posterior approximation randomly proposes a (negated) final momentum from the proposal $r_v$, and runs the leapfrog integrator to find a plausible initial location $x_0$. Finally, the (deterministic) meta-meta-posterior finds the initial momentum that could have taken $x_0$ to $x$.

§.§ Antithetic Sampling

Consider a target $\tilde{\pi}(x)$ and a proposal $q(x)$ that approximates $\pi$. Suppose $q$ is invariant under some bijective transformation $T$:
$$\forall x, q(x) = q(T(x)).$$
For example, a univariate Gaussian proposal with mean $\mu$ is invariant under $T(x) = 2\mu - x$. Antithetic sampling generates a sample $x$ from $q$, but instead of using the estimator $\hat{Z} = \tilde{\pi}(x)/q(x)$, it uses
$$\hat{Z} = \frac{\tilde{\pi}(x) + \tilde{\pi}(T(x))}{2q(x)}.$$

RAVI Inference Strategy: Antithetic Sampling
antitheticantithetic($\tilde{\pi}, q, T$).qantitheticmantithetic($\tilde{\pi}, q, T$).M($x$).q
infalgPosterior Approx.
InfersTarget of inference
AuxAuxiliary variables

latent variable $x$
sampled $x_0$, choice $b$
$x_0 \sim q$
$w_0 \gets \tilde{\pi}(x_0) / q(x_0)$
$w_1 \gets \tilde{\pi}(T(x_0)) / q(x_0)$
$b \sim \texttt{Bernoulli}(\frac{w_1}{w_0 + w_1})$
$bT(x_0) + (1-b)x_0$

metaalgMeta-Posterior Approx.

sampled $x_0$, choice $b$
$b \sim \text{Bernoulli}(0.5)$
$x_0 \gets bT(x) + (1-b)x$
$(x_0, b)$

This can be justified as Algorithm 1 () applied to the strategy [infstrat:antithetic]. The posterior approximation generates an initial sample $x_0 \sim q$, evaluates both $x_0$ and $T(x_0)$ as possible proposals, and selects one. The meta-posterior approximation must recover whether $x$ or its transformed version was the sampled one; it does so by flipping a fair coin, which is optimal when $T = T^{-1}$, i.e., when $T$ is an involution. In the general case a lower-variance estimator could be derived by setting $\mathcal{M}(x).q$ to the exact posterior of the proposal process. Antithetic sampling can also be generalized to the case where a finite family of bijective transformations $T_i$ are available.

Note that although the final expression for $\hat{Z}$ falls out of this inference strategy only when $q(x) = q(T(x))$ for all $x$, nothing in the inference strategy itself exploits this assumption, and the same inference strategy could be applied to $T$ without this property, to derive other estimators that—intuitively—simultaneously consider a proposal $x$ and a deterministic function of it $T(x)$ as possible locations.

§ ABSOLUTE CONTINUITY

When we defined inference strategies $\mathcal{S}$ targeting $\pi$, we required that $\mathcal{S}.q$ and $\pi$ be mutually absolutely continuous, a stronger requirement than in importance sampling. We now consider relaxing this requirement, by requiring only one-sided absolute continuity. We define two kinds of inference strategy, depending on which direction of absolute continuity holds:

* An inference strategy $\mathcal{S}$ targeting $\pi$ is wide if $\pi$ is absolutely continuous with respect to $\mathcal{S}.q$, and either $\mathcal{S}.q$ has a tractable marginal density or $\mathcal{S}.\mathcal{M}(x)$ is a narrow inference strategy targeting $\mathcal{S}.q(\cdot \mid x)$ for all $x$.

* An inference strategy $\mathcal{S}$ targeting $\pi$ is narrow if $\mathcal{S}.q$ is absolutely continuous with respect to $\pi$, and either $\mathcal{S}.q$ has a tractable marginal density or $\mathcal{S}.\mathcal{M}(x)$ is a wide inference strategy targeting $\mathcal{S}.q(\cdot \mid x)$ for all $x$.

Then an inference strategy as defined in the main paper is one that is both wide and narrow.

Narrow inference strategies can serve as variational families within variational inference algorithms. Wide inference strategies can be used as importance sampling and SMC proposals, as well as variational families for amortized variational inference.  Inference strategies used as MCMC proposals must be both wide and narrow.

§ OTHER APPLICATIONS OF RAVI INFERENCE STRATEGIES

§.§ Rejection sampling with RAVI

As in any properly weighted sampler, if the weights produced by Alg. 1 can be bounded above by a constant $M$, a RAVI inference strategy can be used for exact inference via rejection sampling: a sample $(x, \hat{Z})$ is drawn using Alg. 1, and then accepted with probability $\frac{\hat{Z}}{M}$. The weight $\hat{Z}$ for an inference strategy can be viewed as a product of the normalizing constant $Z$ with normalized importance weights $w_\mathcal{S} = \frac{\pi(x)}{\mathcal{S}.q(x)}$, $w_{\mathcal{S}.\mathcal{M}(x)} = \frac{\mathcal{S}.q(r \mid x)}{\mathcal{S}.\mathcal{M}(x).q(r)}$, and so on. As such, if upper bounds $M_Z$ and $M_\mathcal{S}$, $M_{\mathcal{S}.\mathcal{M}(x)}$, etc. can be found for these quantities, the product of these bounds is a bound on $\hat{Z}$. Thus, as in properly weighted sampling and in variational inference with RAVI, it is possible to reason about the RAVI inference strategy compositionally, in terms of bounds at each layer of nesting.

§.§ Estimating KL divergences between models with RAVI inference strategies equipped

Suppose $p(y) = \int p(x, y) \text{d}x$ and $q(y) = \int q(z, y) \text{d}z$ are mutually absolutely continuous distributions over some space $\mathcal{Y}$. Suppose also that we have two families of inference strategies, $\mathcal{S}_p(y)$ and $\mathcal{S}_q(y)$, targeting $p(x\mid y)$ and $q(z \mid y)$ respectively. Then the AIDE algortihm [Cusumano-Towner and Mansinghka, 2017] can be adapted to give a stochastic upper bound on the symmetric KL divergence between $p(y)$ and $q(y)$.

First, we generate $(x, y_p) \sim p$, $(z, y_q) \sim q$, and run $\texttt{HME}$ on each pair to obtain weights $w^p_p$ and $w^q_q$ respectively. Then, we run $\texttt{IMPORTANCE}$ on $p$ with data $y_q$, and on $q$ with data $y_p$, to obtain weights $w^p_q$ and $w^q_p$ respectively. Finally, we sum the logs of the foru weights, to give an estimate $\hat{D}$ whose expectation is:

\mathbb{E}[\hat{D}] =
\mathbb{E}_{y \sim p}[\UU(p, y, \mathcal{S}_p(y)) - \LL(q, y, \mathcal{S}_q(y))] + \mathbb{E}_{y \sim q}[\UU(q, y, \mathcal{S}_q(y)) - \LL_p(p, y, \mathcal{S}_p(y))] \geq KL(p || q) + KL(q || p).

As the marginal likelihood bounds $\UU$ and $\LL$ become tighter, this expectation approaches the true symmetric KL between $p$ and $q$, i.e., $D = KL(p|| q) + KL(q || p)$. Theorem 4 allows us to characterize the tightness of these bounds, and thus of the stochastic upper bound $\hat{D}$ on the symmetric KL, in terms of KL divergences between successive layers of each inference strategy. Improving inference at any layer of the inference strategy tightens the bound $\hat{D}$, yielding less biased estimates of $D$.

§ REPARAMETERIZATION TRICK GRADIENT ESTIMATORS

In this section, we present versions of Algorithms 3 and 4 that
utilize reparameterization gradients, rather than score function
gradients. Using these algorithms requires that an inference
strategy be reparameterizable.

Definition: A reparameterizable inference strategy $\mathcal{S}$ with arguments $\theta$ specifies:

* A reparameterizable posterior approximation $\mathcal{S}.q$, which is one of:

* a tractable proposal: a tuple $(\mathcal{S}.q(x; \theta), \mathcal{S}.q.g(\epsilon), \mathcal{S}.q.f(\epsilon, \theta)$, such that $q$ is the pushforward of $g$ by $f$; or

* an intractable proposal: a tuple $(\mathcal{S}.q(r, x; \theta), \mathcal{S}.q.g(\epsilon_r, \epsilon_x), \mathcal{S}.q.f_r(\epsilon_r, \theta), \mathcal{S}.q.f_x(\epsilon_x, \theta))$, such that $q$ is the pushforward of $g$ by $\lambda (\epsilon_r, \epsilon_x). (f_r(\epsilon_r, \theta), f_x(\epsilon_x, \theta))$.

* If the latter, a reparameterizable meta-inference strategy $\mathcal{S}.\mathcal{M}$, with arguments $(x, \theta)$, that given argument $(x, \theta)$, targets $\mathcal{S}.q(r \mid x; \theta)$.

Now, reparameterized estimators can be derived by applying standard automatic differentiation to the following algorithm, which only samples from distributions that do not depend on parameters:

Algorithm 6: RAVI ELBO estimator ($\texttt{ELBO}$)
unnormalized model $\tilde{p}(x)$
inference strategy $\mathcal{S}$ with arguments
arguments $\theta$
unbiased estimates of $\mathcal{L}$ (differentiable w.r.t. $\theta$)
$\mathcal{S}.q$ has a tractable marginal density
$\epsilon_x \sim \mathcal{S}.q.g$
$x \gets \mathcal{S}.q.f(\epsilon_x, \theta)$
$\hat{U} \gets \log \mathcal{S}.q(x; \theta)$

$\mathcal{S}.q(x; \theta) = \int \mathcal{S}.q(r, x; \theta)\text{d}r$
$(\epsilon_r, \epsilon_x) \sim \mathcal{S}.q.g$
$(x, r) \gets (\mathcal{S}.q.f_x(\epsilon_x, \theta), \mathcal{S}.q.f_r(\epsilon_r, \theta))$
$\hat{U} \gets \texttt{EUBO}(\mathcal{S}.q(\cdot, x; \theta), r, \mathcal{S}.\mathcal{M}, (x, \theta))$

$\log \tilde{p}(x) - \hat{U}$

Algorithm 7: RAVI EUBO estimator ($\texttt{EUBO}$)
unnormalized model $\tilde{p}(x)$
exact sample $x \sim p(x)$
inference strategy $\mathcal{S}$ with arguments
arguments $\theta$
unbiased estimate of $\mathcal{U}$ (differentiable w.r.t. $\theta$)
$\mathcal{S}.q$ has a tractable marginal density
$\hat{L} \gets \log \mathcal{S}.q(x; \theta)$

$\mathcal{S}.q(x; \theta) = \int \mathcal{S}.q(r, x; \theta)\text{d}r$
$\hat{L} \gets \texttt{ELBO}(\mathcal{S}.q(\cdot, x; \theta), \mathcal{S}.\mathcal{M}, (x, \theta))$

$\log \tilde{p}(x) - \hat{L}$

Note that in fact only every other posterior approximation in the unrolled strategy requires a reparameterized version: Algorithm 7 never samples from its $\mathcal{S}.q$, only evaluates the densities.

It would be interesting to develop variants of these algorithms that allow users to combine score-function and reparameterization estimation at different layers of nesting, or exploit other variance reduction tactics compositionally.