# Interpolating Between Sampling and Variational Inference with Infinite
Stochastic Mixtures

Richard D. Lange Dept. of Neurobiology
University of Pennsylvania
Philadelphia, Pennsylvania, USA  Ari S. Benjamin  Dept. of Neurobiology
University of Pennsylvania
Philadelphia, Pennsylvania, USA  Ralf M. Haefner∗ Dept. of Brain and Cognitive
Sciences
University of Rochester
Rochester, New York, USA  Xaq Pitkow∗ Baylor College of Medicine
Rice University Houston, Texas, USA

###### Abstract

Sampling and Variational Inference (VI) are two large families of methods for
approximate inference that have complementary strengths. Sampling methods
excel at approximating arbitrary probability distributions, but can be
inefficient. VI methods are efficient, but may misrepresent the true
distribution. Here, we develop a general framework where approximations are
stochastic mixtures of simple component distributions. Both sampling and VI
can be seen as special cases: in sampling, each mixture component is a delta-
function and is chosen stochastically, while in standard VI a single component
is chosen to minimize divergence. We derive a practical method that
interpolates between sampling and VI by solving an optimization problem over a
mixing distribution. Intermediate inference methods then arise by varying a
single parameter. Our method provably improves on sampling (reducing variance)
and on VI (reducing bias+variance despite increasing variance). We demonstrate
our method’s bias/variance trade-off in practice on reference problems, and we
compare outcomes to commonly used sampling and VI methods. This work takes a
step towards a highly flexible yet simple family of inference methods that
combines the complementary strengths of sampling and VI.

##  1 Introduction

Figure 1: Conceptual introduction on a toy 2D example. a) Sampling methods
approximate the underlying ${\rm p}(\mathbf{x})$ with a stochastic set of
representative points. b) Variational Inference (VI) methods begin by
selecting an approximating distribution family, ${\rm q}(\mathbf{x};\theta)$,
here an isotropic Gaussian plotted as an ellipse at its $1\sigma$ contour. The
optimal parameters $\theta^{*}$ are chosen to minimize ${\rm KL}({\rm
q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))$. c) We propose using a stochastic
mixture of component distributions, where _parameters_ $\theta$ are sampled
rather than the variable of interest $\mathbf{x}$.

We are concerned with the familiar and general case of approximating a
probability distribution, such as occurs in Bayesian inference when both the
prior over latent variables and the likelihood function connecting them to
data are known, but computing the posterior exactly is intractable. There are
two largely separate families of techniques for approximating such intractable
inference problems: Markov Chain Monte Carlo (MCMC) sampling, and Variational
Inference (VI) [Bishop, 2006, Murphy, 2012].

Sampling-based methods, including MCMC, approximate a distribution with a
finite set of representative points. MCMC methods are stochastic and
sequential, generating a sequence of sample points that, given enough time,
become representative of the underlying distribution increasingly well. MCMC
sampling is (typically) asymptotically unbiased, at the expense of high
variance, leading to long run times in practice. Similar to the approach we
take here, sampling methods are studied at different scales: both in terms of
their asymptotic limit (i.e. their bias at infinitely many samples) and their
practical behavior for finite samples or other resource limits [Korattikara et
al., 2014, Angelino et al., 2016].

Variational Inference (VI) refers to methods that produce an approximate
distribution by minimizing some quantification of divergence between the
approximation and the desired posterior distribution [Blei et al., 2017, Zhang
et al., 2019]. For the purposes of this paper, we will use VI to refer to the
most common flavor of variational methods, namely minimizing the Kullback-
Leibler (${\rm KL}$) divergence between an approximate distribution from a
fixed family and the desired distribution [Bishop, 2006, Wainwright and
Jordan, 2008, Murphy, 2012, Blei et al., 2017]. The best-fitting approximate
distribution is often used directly as a proxy for the true posterior in
subsequent calculations, which can greatly simplify those downstream
calculations if the approximate distribution is itself easy to integrate. In
contrast to MCMC, VI is often used in cases where speed is more important than
asymptotic bias [Angelino et al., 2016, Blei et al., 2017, Zhang et al.,
2019].

In this work, our goal is to develop an intermediate family of methods that
“interpolate” between MCMC and VI, inspired by a simple and intuitive picture
(Figure 1): we propose applying sampling methods _in the space of variational
parameters_ such that the resulting approximation is a stochastic mixture of
variational “component” distributions [Yin and Zhou, 2018]. This extends
sampling by replacing the sampled points with extended components, and it
extends VI by replacing the single best-fitting variational distribution with
a stochastic mixture of more localized components. This is qualitatively
distinct from previous variational methods that use _stochastic optimization_
: rather than stochastically optimizing a single variational approximation
[Hoffman et al., 2013, Salimans et al., 2015], we use stochasticity to
construct a _random mixture_ of variational components that achieves lower
asymptotic bias than any one component could. As we will show below, this
framework generalizes both sampling and VI, where sampling and VI emerge as
special cases of a single optimization problem.

This paper is organized as follows. In section 2, we set up the problem and
our notation, and describe how both classic sampling and classic VI can be
understood as special cases of stochastic mixtures. In section 3, we introduce
an intuitive framework for reasoning about infinite stochastic mixtures and
define an optimization problem that captures the trade-off between sampling
and VI. Section 4 introduces an approximate objective and closed-form solution
and describes a simple practical algorithm. Section 5 gives empirical and
theoretical results that show how our method interpolates the bias and
variance of sampling and VI. Finally, section 6 concludes with a summary,
related work, limitations, and future directions.

##  2 SETUP AND NOTATION

Let ${\rm p}^{*}(\mathbf{x})={\rm Z}{\rm p}(\mathbf{x})$ denote the
unnormalized probability distribution of interest, with unknown normalizing
constant ${\rm Z}$. For instance, in the common case of a probabilistic model
with latent variables $\mathbf{x}$, observed data $\mathcal{D}$, and joint
distribution ${\rm p}(\mathbf{x},\mathcal{D})$, we are interested in
approximations to the posterior distribution ${\rm
p}(\mathbf{x}|\mathcal{D})$. This is intractable in general, but we assume
that we have access to the un-normalized posterior ${\rm
p}^{*}(\mathbf{x}|\mathcal{D})=\frac{1}{{\rm Z}}{\rm
p}(\mathcal{D}|\mathbf{x}){\rm p}(\mathbf{x})$.111To reduce clutter,
$\mathcal{D}$ will be dropped in the remainder of the paper, and we will use
only ${\rm p}(\mathbf{x})$ and ${\rm p}^{*}(\mathbf{x})$. Let ${\rm
q}(\mathbf{x};\theta)$ be any “simple” distribution that may be used used in a
classic VI context (such as mean-field or Gaussian), and let ${\rm
m}_{T}(\mathbf{x})$ be a mixture containing $T$ of these simple distributions
as components, defined by a set of $T$ parameters
$\\{{\theta^{(1)},\ldots,\theta^{(T)}\\}}$:

${\rm m}_{T}(\mathbf{x})\equiv\frac{1}{T}\sum\limits_{t=1}^{T}{\rm
q}(\mathbf{x};\theta^{(t)})\,.$ (1)

For example, if ${\rm q}$ is a multivariate normal with mean $\mu$ and
covariance $\Sigma$, then $\theta^{(t)}=\\{{\mu^{(t)},\Sigma^{(t)}}\\}$ and
${\rm m}_{T}(\mathbf{x})$ would be a mixture of $T$ component normal
distributions [Gershman et al., 2012].

We will study properties of distributions over component parameters, which we
denote $\psi(\theta)$ [Ranganath et al., 2016]. If the set of $\theta^{(t)}$
is drawn randomly from $\psi(\theta)$, then as $T\rightarrow\infty$, ${\rm
m}_{T}(\mathbf{x})$ approaches the idealized infinite mixture,

${\rm m}(\mathbf{x})\equiv\int_{\theta}{\rm
q}(\mathbf{x};\theta)\psi(\theta){\rm d}\theta\,.$ (2)

#### Sampling and VI as special cases of the mixing distribution.

Let $\theta^{*}=\operatorname*{arg\,min}_{\theta}{\rm KL}({\rm
q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))$ be the parameters corresponding
to the classic single-component variational solution. VI corresponds to the
special case where the mixing distribution $\psi(\theta)$ is a Dirac delta
around $\theta^{*}$, or $\psi(\theta)=\delta(\theta-\theta^{*})$, in which
case the mixture ${\rm m}_{T}(\mathbf{x})$ is equivalent to ${\rm
q}(\mathbf{x};\theta^{*})$ regardless of the number of components $T$.
Sampling can also be seen as a special case of $\psi(\theta)$ in which each
component narrows to a Dirac delta ($\psi(\theta)$ places negligible mass on
regions of $\theta$-space where components have appreciable width), and the
means of the components are distributed according to ${\rm p}(\mathbf{x})$.
This requires that the component family ${\rm q}(\mathbf{x};\theta)$ is
capable of expressing a Dirac-delta at any point $\mathbf{x}$, such as a
location-scale family. Thus, both sampling and VI can be seen as limiting
cases of stochastic mixture distributions, ${\rm m}_{T}(\mathbf{x})$, defined
by a distribution over component parameters, $\psi(\theta)$. In what follows,
we will show how designing the mixing distribution $\psi(\theta)$ allows us to
create mixtures that trade-off the complementary strengths of sampling and VI.

Figure 2: Left: Understanding mixtures in terms of Mutual Information and
Expected KL. a) The quality of any infinite mixture (in terms of ${\rm
KL}({\rm m}||{\rm p})$) is given by its distance from the y=x line (black
diagonal line). b) Two unreachable regions are shaded in gray: above the y=x
line (because ${\rm KL}({\rm m}||{\rm p})\geq 0$), and to the left of the
single-component variational solution, since VI achieves the minimum ${\rm
KL}({\rm q}||{\rm p})$. c) When $\psi(\theta)=\delta(\theta-\theta^{*})$ as in
classic VI, Expected KL is at its minimum and Mutual Information is zero.
Increasing the expressiveness of ${\rm q}$ corresponds to moving left along
the x-axis (blue arrow). d) Because sampling is unbiased, it is a mixture that
lives on the ${\rm KL}({\rm m}||{\rm p})=0$ or $y=x$ line. If $\mathbf{x}$ is
discrete, the coordinates of the point marked (d) are
$(\mathcal{H}[\mathbf{x}],\mathcal{H}[\mathbf{x}])$, i.e. the entropy of ${\rm
p}(\mathbf{x})$. When $\mathbf{x}$ is continuous, both Mutual Information and
Expected KL grow unboundedly together as the individual components narrow. e)
Any point on the y=x line implies ${\rm m}(\mathbf{x})={\rm p}(\mathbf{x})$,
and this may be possible without resorting to sampling for certain
combinations of ${\rm p}$ and ${\rm q}$. However, such mixtures are not
guaranteed to exist for all problems, and are difficult to find due to the
intractability of Mutual Information. f) We propose a family of mixture
approximations, parameterized by $\lambda$, that connects VI to sampling in a
natural and principled way. Points on this curve correspond to solutions to
the (approximate version of the) objective in (5). Middle: Examples in a 1D
toy problem, where ${\rm p}(\mathbf{x})$ is an unequal mixture of two heavy-
tailed distributions (black lines), and ${\rm q}(\mathbf{x};\theta)$ is a
single Gaussian component with parameters $\theta=\\{{\mu,\log\sigma}\\}$
(transluscent red components). Right: Varying $\lambda$ controls the mixing
distribution over $\theta$ (contours). Red points correspond to the Gaussian
components in the middle.

##  3 Conceptual Framework

###  3.1 Decomposing ${\rm KL}({\rm m}||{\rm p})$ Into Mutual Information and
Expected KL

The idealized infinite mixture ${\rm m}(\mathbf{x})$ is fully defined by the
chosen component family ${\rm q}(\mathbf{x};\theta)$ and the mixing
distribution $\psi(\theta)$. Consider the variational objective with respect
to the entire mixture, ${\rm KL}({\rm m}||{\rm p})$:

$\begin{split}{\rm KL}({\rm m}||{\rm p})=\int_{\mathbf{x}}{\rm
m}(\mathbf{x})\log\frac{{\rm m}(\mathbf{x})}{{\rm p}^{*}(\mathbf{x})}{\rm
d}\mathbf{x}+\log{\rm Z}\,,\end{split}$ (3)

where ${\rm Z}$ is the normalizing constant of ${\rm p}^{*}(\mathbf{x})$ and
is irrelevant for constructing ${\rm m}(\mathbf{x})$. Instead of (3), one can
use the equivalent objective of maximizing the _Evidence Lower BOund_ or ELBO
[Bishop, 2006, Murphy, 2012, Blei et al., 2017]. Regardless, minimizing (3) or
maximizing the ELBO for mixtures is intractable in general. However, as first
shown by Jaakkola and Jordan [1998] for finite mixtures, it admits the
following useful decomposition:

$\begin{split}{\rm KL}({\rm m}||{\rm
p})&=\underbrace{\int_{\theta}\psi(\theta)\int_{\mathbf{x}}{\rm
q}(\mathbf{x};\theta)\log\frac{{\rm q}(\mathbf{x};\theta)}{{\rm
p}^{*}(\mathbf{x})}{\rm d}\mathbf{x}\,{\rm d}\theta}_{\text{(i) Expected
KL}}\\\ &-\underbrace{\int_{\theta}\psi(\theta)\int_{\mathbf{x}}{\rm
q}(\mathbf{x};\theta)\log\frac{{\rm q}(\mathbf{x};\theta)}{{\rm
m}(\mathbf{x})}{\rm d}\mathbf{x}\,{\rm d}\theta}_{\text{(ii) Mutual
Information }\mathcal{I}[\mathbf{x};\theta]}\end{split}$ (4)

(dropping $\log{\rm Z}$). The first term, (i), is the Expected KL Divergence
for each component when the parameters are drawn from $\psi(\theta)$. This
term quantifies, on average, how well the mixture components match the target
distribution. In isolation, Expected KL is minimized when all components
individually minimize ${\rm KL}({\rm q}||{\rm p})$, i.e. when
$\psi(\theta)\rightarrow\delta(\theta-\theta^{*})$. This tendency to
concentrate $\psi(\theta)$ to the single best variational solution is balanced
by the second term, (ii), which is the Mutual Information between $\mathbf{x}$
and $\theta$, which we will write $\mathcal{I}[\mathbf{x};\theta]$, under the
joint distribution ${\rm q}(\mathbf{x};\theta)\psi(\theta)$. This term should
be _maximized_ , and, importantly, it does not depend on ${\rm
p}^{*}(\mathbf{x})$. Mutual Information is maximized when the components are
as diverse as possible, which encourages the components to become narrow and
to spread out over diverse regions of $\mathbf{x}$ _regardless_ of how well
they agree with ${\rm p}(\mathbf{x})$. This decomposition of ${\rm KL}({\rm
m}||{\rm p})$ into Mutual Information (between $\mathbf{x}$ and $\theta$) and
Expected KL (between ${\rm q}$ and ${\rm p}$) is convenient because
approximations to Mutual Information are well-studied, and minimizing Expected
KL can leverage standard tools from VI.

###  3.2 Trading Off Between Mututal Information and Expected KL

We will refer back to this decomposition of the ${\rm KL}({\rm m}||{\rm p})$
objective into Expected KL (between ${\rm q}$ and ${\rm p}$) and Mutual
Information (between each $\mathbf{x}$ and $\theta$) throughout. Figure 2
depicts a two-dimensional space with Expected KL on the x-axis and Mutual
Information on the y-axis. Any given mixing distribution $\psi$ can be placed
as a point in this space, but in general many $\psi$’s may map to the same
point.

Sampling and VI live at extreme points in this space. Classic VI, where
$\psi(\theta)=\delta(\theta-\theta^{*})$, corresponds to the blue point (c),
because by definition $\theta^{*}$ achieves the minimum possible ${\rm KL}$,
and $\mathcal{I}[\mathbf{x};\theta]$ is zero. Classic sampling corresponds to
the green point (d), with $\psi(\theta)$ placing mass only on Dirac-delta-like
components, and selecting each component with probability ${\rm p}(\mu)$,
where $\mu$ is the mean of ${\rm q}$ determined by $\theta$.

Towards the goal of constructing mixtures that trade-off properties of
sampling and VI, we propose to view the two terms in (4) as separate
objectives that may be differently weighted, and maximizing the objective

$\mathcal{L}(\psi,\lambda)=\mathcal{I}[\mathbf{x};\theta]-\lambda\mathbb{E}_{\psi}\left[{\rm
KL}({\rm q}||{\rm p})\right]$ (5)

for a given hyperparameter $\lambda$ with respect to the mixing distribution
$\psi$. This objective may alternatively be viewed as the Lagrangian of a
constrained optimization problem over the mixing density $\psi$, where Mutual
Information is maximized subject to a constraint on Expected KL. This is a
concave maximization problem with linear constraints, defining a Pareto front
of solutions that each achieve a different balance between Expected KL and
Mutual Information. In practice, maximizing Mutual Information necessitates
approximations [Poole et al., 2019], so there may be good mixture
approximations that are not found in practice, such as the yellow point (e) in
Figure 2. In section 4 below, we use an approximation to Mutual Information
that has the property, illustrated by the orange curve (f) in Figure 2, of
connecting VI (c) to sampling (d), controlled by varying $\lambda$. As shown
on the right of Figure 2, our method produces mixtures that behave like
classic samples when $\lambda=1$, that behave like classic VI when
$\lambda\rightarrow\infty$, and that exhibit intermediate behavior at
intermediate values of $\lambda$.

We emphasize that this frame is quite general: any stochastic mixture can be
reasoned about in terms of its Expected KL and Mutual Information, and this is
a natural space in which to think about interpolating sampling and VI. A
similar decomposition of ${\rm KL}({\rm m}||{\rm p})$ (or the ELBO) has been
used by previous methods that optimize mixtures [Zobay, 2014, Jaakkola and
Jordan, 1998, Gershman et al., 2012, Yin and Zhou, 2018]. The primary
difference between these previous methods is how they approximate (or lower-
bound) Mutual Information. In the next section, we introduce a new
approximation that is particularly efficient, and is the first to our
knowledge that can produce sampling-like behavior with finitely many
components.

##  4 Approximate objective

Maximizing Mutual Information, as is required by (5), is a notoriously
difficult problem that arises in many domains, and there is a large collection
of approximations and bounds in the literature [Jaakkola and Jordan, 1998,
Brunel and Nadal, 1998, Gershman et al., 2012, Wei and Stocker, 2016,
Kolchinsky and Tracey, 2017, Poole et al., 2019]. Previous work has optimized
_finite_ mixtures by considering how each of $T$ components interacts with the
other $T-1$ components, resulting in quadratic scaling with $T$ [Gershman et
al., 2012, Guo et al., 2016, Miller et al., 2017, Kolchinsky and Tracey, 2017,
Yin and Zhou, 2018, Poole et al., 2019]. Beginning instead with _infinite_
mixtures, we find that the local geometry of $\theta$-space is sufficient to
provide an approximation to Mutual Information _that can be evaluated
independently for each value of $\theta$_.

###  4.1 Stam’s inequality

Mutual Information between $\mathbf{x}$ and $\theta$ can be written as

$\displaystyle\mathcal{I}[\mathbf{x};\theta]$
$\displaystyle=\mathcal{H}[\theta]-\mathbb{E}_{{\rm
m}(\mathbf{x})}\left[\mathcal{H}[\hat{\theta}|\mathbf{x}]\right]$
$\displaystyle=\mathcal{H}[\theta]-\mathbb{E}_{\psi(\theta)}\big{[}\underbrace{\mathbb{E}_{{\rm
q}(\mathbf{x}|\theta)}[\mathcal{H}[\hat{\theta}|\mathbf{x}]]}_{\mathcal{H}[\hat{\theta}|\theta]}\big{]}$
(6)

where $\mathcal{H}[\theta]$ is the entropy of $\psi(\theta)$ and
$\mathcal{H}[\hat{\theta}|\mathbf{x}]$ is the entropy of ${\rm
q}(\hat{\theta}|\mathbf{x})=\frac{{\rm
q}(\mathbf{x};\hat{\theta})\psi(\hat{\theta})}{{\rm m}(\mathbf{x})}$, i.e. the
distribution of _inferred_ $\theta$ values for a given $\mathbf{x}$. The
second line follows simply from expanding the definition of ${\rm
m}(\mathbf{x})$ in the outer expectation. The term
$\mathcal{H}[\hat{\theta}|\theta]$ can be thought of in terms of a statistical
estimation problem: $\hat{\theta}$ is the “recovered” value of $\theta$ after
passing through the “channel” $\mathbf{x}$. Bounding the error of such
estimators is a well-studied problem in statistics.

From (6), a lower-bound on Mutual Information can be derived from an _upper
bound_ on $\mathcal{H}[\hat{\theta}|\theta]$ for each $\theta$. For this, we
draw inspiration from Stam’s inequality [Stam, 1959, Dembo et al., 1991, Wei
and Stocker, 2016], which states

$\mathcal{H}[\hat{\theta}|\theta]\leq\frac{1}{2}\log\left|2\pi
e\mathcal{F}(\theta)^{-1}\right|\,,$ (7)

where $|\cdot|$ is a determinant, and $\mathcal{F}(\theta)$ is the Fisher
Information Matrix, defined as

$\displaystyle\mathcal{F}(\theta)_{ij}=-\mathbb{E}_{{\rm
q}(\mathbf{x};\theta)}\left[\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\log{\rm
q}(\mathbf{x};\theta)\right]\,.$

The Fisher Information Matrix is also the local metric on the _statistical
manifold_ with coordinates $\theta$ [Amari, 2016]; it is used to quantify how
“distinguishable” $\theta$ is from $\theta+d\theta$. Note that (7) can be
viewed as the entropy of a Gaussian approximation to ${\rm
q}(\hat{\theta}|\mathbf{x})$ with precision matrix $\mathcal{F}(\theta)$; this
approximation is most accurate when ${\rm q}(\mathbf{x};\theta)$ itself is
narrow and approximately Gaussian [Wei and Stocker, 2016].

Combining (6) and (7), we propose to use

$\mathcal{I}_{\mathcal{F}}[\mathbf{x};\theta]\equiv\mathcal{H}[\theta]-\frac{1}{2}\mathbb{E}_{\psi(\theta)}\left[\log\left|2\pi
e\mathcal{F}(\theta)^{-1}\right|\right]$ (8)

as a proxy for the intractable $\mathcal{I}[\mathbf{x};\theta]$ in (5).

Note that $\mathcal{I}_{\mathcal{F}}[\mathbf{x};\theta]$ is not strictly a
_bound_ on $\mathcal{I}[\mathbf{x};\theta]$, but may be seen as an
_approximation_ to it [Wei and Stocker, 2016]. Briefly, this is because the
original Stam’s inequality, as stated in (7), assumes $\theta$ is a scalar
location parameter, and assumes the high-precision limit where ${\rm
q}(\hat{\theta}|\mathbf{x})$ is well-approximated by a Gaussian. Despite this,
$\mathcal{I}_{\mathcal{F}}[\mathbf{x};\theta]$ is well-suited for our
purposes, since (i) it leads to a remarkably simple and easy to implement
expression for $\psi(\theta)$ below; (ii) we can prove that it leads to
sampling when $\lambda=1$ and VI when $\lambda\rightarrow\infty$; and (iii)
the inequality in (7) is nonetheless likely to be strict, since we neglect the
prior information contained in $\psi(\theta)$ when estimating $\hat{\theta}$
and therefore over-estimate the entropy.222By analogy to the Bayesian Cramér-
Rao bound [Gill and Levit, 1995, Fauß et al., 2021], a tighter variant of (7)
could be derived that takes into account the prior, though possibly at the
expense of added complexity; we leave this to future work.

###  4.2 Closed-form mixing distribution

Substituting $\mathcal{I}_{\mathcal{F}}[\mathbf{x};\theta]$ for
$\mathcal{I}[\mathbf{x};\theta]$ in (5) gives the following approximate
objective,

$\mathcal{L}_{\mathcal{F}}(\psi,\lambda)=\mathcal{H}[\theta]+\mathbb{E}_{\psi}\left[\frac{1}{2}\log|\mathcal{F}|-\lambda\,{\rm
KL}({\rm q}||{\rm p}^{*})\right]$ (9)

having dropped additive constants and using
$\log|\mathcal{F}^{-1}|=-\log|\mathcal{F}|$. This now resembles a maximum-
entropy problem with an expected-value constraint, which has the following
simple closed-form solution:

$\log\psi(\theta)=\frac{1}{2}\log|\mathcal{F}(\theta)|-\lambda\,{\rm KL}({\rm
q}(\mathbf{x};\theta)||{\rm p}^{*}(\mathbf{x}))$ (10)

again dropping additive constants. Equation (10) is strikingly simple, and
amenable to many existing MCMC sampling methods for drawing samples of
$\theta$ from $\psi$.

Despite being derived from an approximation to our original objective, (10)
nonetheless contains both sampling and VI as special cases. As
$\lambda\rightarrow\infty$, the ${\rm KL}$ term dominates and $\psi(\theta)$
concentrates to $\delta(\theta-\theta^{*})$, reproducing VI. When $\lambda=1$,
this mixing distribution also corresponds to “sampling” in the following
sense:

######  Definition 1 (Sampling)

A stochastic mixture, defined by the component family ${\rm
q}(\mathbf{x};\theta)$ and mixing distribution $\psi(\theta)$, is considered
to be “sampling” if (i) it is unbiased in the limit of infinitely many
components, i.e. ${\rm m}(\mathbf{x})\rightarrow{\rm p}(\mathbf{x})$; and,
(ii) it consists of non-overlapping components. That is, for small values of
of $0<\epsilon\ll 1$, wherever ${\rm q}(\mathbf{x};\theta_{i})>\epsilon$, with
high probability, ${\rm q}(\mathbf{x};\theta_{j})<\epsilon$, for all pairs
$\theta_{i},\theta_{j}$ drawn independently from $\psi(\theta)$.

Lemma 4 in Appendix A.2 establishes that $\psi(\theta)$ with $\lambda=1$ leads
to sampling as defined here, assuming mixture components ${\rm q}$ are
Gaussian. However, we conjecture that sampling arises from a broader class of
${\rm q}$ components as well.

###  4.3 Implementation

We implemented (10) in Stan [Carpenter et al., 2017], an open-source framework
for probabilistic models and approximate inference algorithms. We sampled
$\theta$ from $\psi(\theta)$ using Stan’s default implementation of the No
U-Turn Sampler (NUTS) [Hoffman and Gelman, 2014], but we emphasize that
samples can be drawn from (10) using any existing sampling method. All
comparisons to existing methods were with Stan’s built-in NUTS sampler (over
$\mathbf{x}$) and its built-in mean-field VI [Kucukelbir et al., 2017].

##  5 Navigating Bias/Variance Trade-Offs For Finite $T$

Figure 3: $\lambda$ controls a bias/variance tradeoff, interpolating between
sampling and VI. a) For an example 2D distribution (the banana distribution),
we set ${\rm q}$ to Gaussian with diagonal covariance and sampled
$\theta\sim\psi(\theta)$ using NUTS (see Appendix B.2 for sampling details).
b) We selected $f(\mathbf{x})$ as a random mixture of sinusoids at different
frequencies. We then calculated the bias and variance of computing
$\mathbb{E}_{{\rm m}_{T}}[f(\mathbf{x})]$. c) The green point and error bars
(“Sampling”) indicate the estimated value of $\mathbb{E}_{{\rm
p}(\mathbf{x})}[f(\mathbf{x})]$ and its variance using NUTS to draw samples of
$\mathbf{x}$. The blue point (“VI”) shows the value of $\mathbb{E}_{{\rm
q}(\mathbf{x};\theta^{*})}[f(\mathbf{x})]$ using Stan’s built-in VI. Our
method is shown in the middle across a range of $\lambda$ values. Low
$\lambda$ provides unbiased but high variance estimators, while high $\lambda$
provides a bias near that of standard VI and a vanishing variance. In panel
(c), we used $T=10$ independent samples for both classic NUTS and our method.
Figure 4: Flexibly trading bias for variance matters for integrating
functions. a) We generated random functions of a specified smoothness by
varying the decay of its power spectrum while randomizing phase. b) The
$\lambda$ with the smallest MSE for a fixed number samples $T=100$ depends on
the integrand’s smoothness. c) Surface plots of MSE for $T=100$ samples,
varying $\lambda$ and $\alpha$. d) Variance is higher for smaller $\lambda$
and more wiggly integrands. e) Bias vanishes near $\lambda=1$.

###  5.1 Reducing Mean Squared Error (MSE)

In this section, we expound the sense in which our method “interpolates”
sampling and VI in terms of bias and variance, both analytically and
empirically. In our experiments, we quantify bias and variance in terms of the
Mean Squared Error (MSE) of the expectation of an arbitrary $f(\mathbf{x})$
using a random mixture of $T$ components, ${\rm m}_{T}(\mathbf{x})$. In Figure
3, we show empirically that by increasing $\lambda$ one can interpolate
between the zero bias but high variance solution, equivalent to sampling, and
the zero variance but high bias solution, equivalent to VI. Between these
extremes, our method smoothly interpolates both bias and variance.

To show this empirically requires choosing a class of functions
$f(\mathbf{x})$. We construct a random smooth function by discrete Fourier
synthesis. Specifically, we select a series of sinusoid plane waves in the
space of $\mathbf{x}$ with increasing frequency $\omega$, random directions
$\mathbf{t}$ and phase $\phi$, such that
$f(\mathbf{x})=\sum_{\omega=1}^{N}a_{\omega}\sin(\omega\mathbf{t}^{\top}\mathbf{x}+\phi_{\omega})$.
The amplitudes $a_{\omega}$ are set according to a power law:
$a_{\omega}=\omega^{-\alpha}$. An example of $f(\mathbf{x})$ is shown in Fig.
3b for $\alpha=-1$, and $\alpha$ is varied in Fig 4. Adjusting $\alpha$ allows
flexibly setting the “wiggliness” of the synthesized function [Stein and
Shakarchi, 2011].

We also tested our algorithm on three reference problems from posteriordb
[Magnusson et al., 2021], now evaluating a large set of random $f$s, defined
on the space of each model’s unconstrained parameters, with $\alpha=-1$
(Figure B.1). The conclusion is similar: across many random $f$s, our
algorithm performs on average as well as or better than both sampling (by
reducing variance) and VI (by reducing bias).

###  5.2 Considerations for selecting $\lambda$

A first practical consideration for the choice of $\lambda$ is the particular
function $f(\mathbf{x})$ to be integrated. Since MSE can be decomposed into
the sum of squared bias and variance, the value of $\lambda$ that minimizes
MSE occurs when
$\frac{\partial\textrm{Bias}^{2}}{\partial\lambda}=-\frac{\partial\textrm{Var}}{\partial\lambda}$.
Any factor that increases the variance but not the bias of an estimate for a
fixed number of components $T$ will push the optimal $\lambda$ towards higher
values.

One such factor is the smoothness of $f(\mathbf{x})$. Classic sampling can
have problematically high variance when $f(\mathbf{x})$ is very jagged, as
single points are not very representative of the surrounding function.
Intuitively, then, we should expect that higher $\lambda$ (more VI-like
mixtures) is preferred when $f(\mathbf{x})$ is more “wiggly.” To show this, we
generated a random function with varying smoothness, integrated it over random
mixtures approximating the 2D banana distribution, and plotted the resulting
MSE, bias, and variance (Fig. 4). We adjusted smoothness by varying the power
law decay, $\alpha$, for a fixed set of phases and wave directions. At any
value of $\lambda$, variance can be seen to increase as $f$ is made more
wiggly. With all else held equal, it is better to trade some variance for bias
when the integrand changes quickly with $\mathbf{x}$.

Another factor that affects the optimal $\lambda$ is the computational budget.
If time allows a large number of $T$ to be sampled, the optimal $\lambda$ will
approach $1$ with a speed that depends on the particular problem
(specifically, on $\frac{\partial\textrm{Bias}^{2}}{\partial\lambda}$). In our
experiments we set a fixed $T$ to demonstrate our algorithm’s properties.
However, if the number of components is not known in advance, a practitioner
may also decrease $\lambda$ adaptively over time as sampling continues.

###  5.3 Analytical results

While the MSE of the expected value of some $f(\mathbf{x})$ is a useful way to
compare approximate inference methods, it depends on the somewhat arbitrary
choice of $f$, and in practice, the $f$’s of interest are often not known at
the time of inference. This motivates using the following alternative
definition of error that is independent of $f$ and closely related to the
variational objective of minimizing ${\rm KL}$ divergence:

$\begin{split}\text{KL error}=\mathbb{E}[{\rm KL}({\rm
m}_{T}(\mathbf{x})||{\rm p}(\mathbf{x}))]=\\\ \underbrace{{\rm KL}({\rm
m}(\mathbf{x})||{\rm p}(\mathbf{x}))}_{\text{KL
bias}}+\underbrace{\mathbb{E}\left[{\rm KL}({\rm m}_{T}(\mathbf{x})||{\rm
m}(\mathbf{x}))\right]}_{\text{KL variance}}\,.\end{split}$ (11)

That is, KL bias is the ${\rm KL}$ divergence from the infinite mixture ${\rm
m}(\mathbf{x})$ to the true distribution, and KL variance is the average ${\rm
KL}$, over realizations of $T$ independent mixture components, from ${\rm
m}_{T}(\mathbf{x})$ to the infinite mixture ${\rm m}(\mathbf{x})$. Note that
KL bias is identical to the infinite-mixture objective we started with in (4).

The following theorem establishes that for all finite $T$, we can always
reduce the KL error, relative to sampling, using some $\lambda>1$.

######  Theorem 1 (Improve on sampling)

If a mixture is sampling as in Definition 1, then $\frac{{\rm d}}{{\rm
d}\lambda}\text{KL bias}=0$ and $\frac{{\rm d}}{{\rm d}\lambda}\text{KL
variance}<0$. Thus, $\frac{{\rm d}}{{\rm d}\lambda}\text{KL error}<0$.

This theorem establishes the intuitive result that the variance of sampling
can be reduced, minimally impacting its bias, by replacing samples with narrow
mixture components. Importantly, Theorem 1 is based on how $\psi(\theta)$
changes with $\lambda$ when using the closed-form expression for $\psi$ we
derived based on the approximate $\mathcal{L}_{\mathcal{F}}$ objective.
Because the theorem is phrased in conditional terms (“if the mixture is
sampling, then…”), we must further show that both conditions of “sampling”
(Definition 1) are met when $\lambda=1$. This is proved in Lemma 4 in Appendix
A.2 for Gaussian components, though we suspect it holds for other component
families as well.

We can also improve on VI using stochastic mixtures. However, this result is
slightly more subtle, as there are three cases where one should expect VI to
be optimal. First, if ${\rm q}$ is in the same family as ${\rm p}$, then ${\rm
q}(\mathbf{x};\theta^{*})={\rm p}(\mathbf{x})$, then is no benefit to
increasing $T$, and reducing $\lambda$ only adds variance. Second, if $T$ is
small – in the most extreme case, if $T=1$ – then reducing $\lambda$ will
again only add variance without reducing bias. Third, if ${\rm p}$ is lighter-
tailed than ${\rm q}(\mathbf{x};\theta^{*})$, then a mixture of nearby ${\rm
q}$s will add variance to ${\rm m}$ [Lindsay, 1983], making the match to ${\rm
p}$ worse. With these three exceptions in mind, the following theorem
establishes conditions where we expect to reduce KL error relative to VI by
using a large but finite $\lambda<\infty$.

######  Theorem 2 (Improve on VI)

Assume that ${\rm p}(\mathbf{x})$ is heaver-tailed than ${\rm
q}(\mathbf{x};\theta^{*})$ and that $\lambda$ is large. Then, there exists
some finite $T_{0}>1$ such that for all $T\geq T_{0}$, $\frac{{\rm d}}{{\rm
d}\lambda}\text{KL error}>0$. Proof: see Appendix A.3.

Note that this result depends on an additional conjecture that relates the
curvature in parameter space of ${\rm KL}({\rm q}||{\rm p})$ to the curvature
of ${\rm KL}({\rm q}||{\rm q}^{*})$ that we believe holds as long as ${\rm p}$
is heavier-tailed than ${\rm q}$. For details, see Appendix A.3.

##  6 DISCUSSION

#### Summary:

Our work provides a new perspective on the relationship between the two
dominant frameworks for approximate inference – sampling and VI – by viewing
both as special cases of inference using a broader class of stochastic
mixtures. Our main theoretical contribution is the framework shown in Figure
2, where mixtures that “interpolate” sampling and VI are analyzed in terms of
how they trade off Mutual Information and Expected KL. We then derived an
easy-to-use method based on an approximation to Mutual Information that uses
the local geometry of the space of variational parameters. To demonstrate the
ease and effectiveness of our method, we implemented it in the popular Stan
language and demonstrated using a small set of reference problems how we
“interpolate” sampling and VI by varying a single parameter, $\lambda$.
Finally, we showed why such an intermediate inference scheme is useful in
practice. On one hand, we proved that it is always possible to improve on
classic sampling ($\lambda=1$) by increasing $\lambda$: our method provably
reduces the variance of sampling while minimally impacting its bias. On the
other hand, our method provably reduces the bias of VI under certain
conditions (and improves overall error if the number of mixture components is
sufficiently large).

#### Time and space complexity:

By approximating Mutual Information using only _local_ geometric information
in (8), in our method each component can be selected independently of the
others. This means we can select and evaluate $T$ components in
$\mathcal{O}(T)$ time and either $\mathcal{O}(T)$ space (if all are stored) or
$\mathcal{O}(1)$ space (if components are evaluated online) – identical to
traditional MCMC sampling algorithms. Further, we can run independent chains
sampling $\theta\sim\psi(\theta)$ for a constant factor speedup. This improves
on past work using mixture approximations, which incurred $\mathcal{O}(T^{2})$
time and $\mathcal{O}(T)$ space complexity, since the optimization problem for
the $T$th component depends on the location of the other $T-1$ components, all
of which must be in memory at once [Jaakkola and Jordan, 1998, Gershman et
al., 2012, Salimans et al., 2015, Guo et al., 2016, Miller et al., 2017,
Acerbi, 2018, Yin and Zhou, 2018] (but the $\mathcal{O}(T^{2})$ complexity may
be hardware-accelerated).

#### Related Work:

The trade-offs between sampling and VI are well-studied, and many methods have
been proposed to “close the gap” between them (see [Angelino et al., 2016,
Zhang et al., 2019] for general reviews). Like these other methods, we aim to
provide good approximations with high computational efficiency and low
variance.

There are many methods that use mixture models to reduce the bias of
variational inference. Theorem 2 shows that our method only “beats” classic VI
when $T>T_{0}$ for some finite but potentially large $T_{0}$. This is the
price we pay for drawing mixture components stochastically [Yin and Zhou,
2018]. When a mixture of $T$ components is _optimized_ rather than _sampled_ ,
bias is reduced and variance remains near zero, as in previous work [Jaakkola
and Jordan, 1998, Gershman et al., 2012, Zobay, 2014, Guo et al., 2016, Miller
et al., 2017], but in previous work this optimization has incurred a
$\mathcal{O}(T^{2})$ cost while our method is $\mathcal{O}(T)$ and can be
further parallelized. Further, with some notable exceptions [Anaya-Izquierdo
and Marriott, 2007, Salimans et al., 2015], most mixture VI methods make
strong assumptions about the family of components [Jaakkola and Jordan, 1998,
Gershman et al., 2012, Acerbi, 2018, Miller et al., 2017]. Our framework and
method is somewhat agnostic to the family of ${\rm q}$, though we have only
rigorously proved that is asymptotically unbiased when using Gaussian
components.

Many methods use sampling in the service of variational inference, or vice
versa, but do not provide a unifying approach to both. These typically use the
samples to compute expectations used to update a variational approximation
[Acerbi, 2018, Miller et al., 2017, Kucukelbir et al., 2017], rather than to
generate the mixture components themselves.

There is also a large number of sampling approaches that aim to improve the
efficiency of sampling by reducing its variance at the cost of some bias. Some
of these use variational approaches as proposal distributions, but ultimately
the posterior is approximated by a set of (possibly weighted) samples of the
latent variables [de Freitas et al., 2001, Korattikara et al., 2014, Ma et
al., 2015, Zhang et al., 2021]. By expanding each sample to a distribution,
our approach allows each sample to cover more space with less variance and
greater efficiency [Nalisnick and Smyth, 2017].

Despite some high-level similarities to other approaches, our framework is
unusual in approximating the posterior by a sampled mixture of variational
approximations. The Mixture Kalman filter [Chen and Liu, 2000] is a special
case of this, which uses a sampled mixture of Gaussians, each constructed as a
Kalman filter. A related approach is to _optimize_ a parameterized function
that generates mixture components [Salimans et al., 2015, Wolf et al., 2016,
Yin and Zhou, 2018], and generative diffusion models can also be seen as a
case of this approach [Sohl-Dickstein et al., 2015, Ho et al., 2020]. Our work
differs in that we derived a closed-form mixing distribution that requires no
additional learning or optimization and that is readily implemented in
existing inference software (Stan, [Carpenter et al., 2017]).

#### Limitations and future work:

Using $\mathcal{I}_{\mathcal{F}}[\mathbf{x};\theta]$ to approximate
$\mathcal{I}[\mathbf{x};\theta]$ reduces the generality of our method, since
the former is most appropriate for narrow and Gaussian-like components [Wei
and Stocker, 2016]. Incorporating prior information from $\psi(\theta)$ into
this bound, generalizing to other kinds of components, or even starting with
alternative bounds on $\mathcal{I}[\mathbf{x};\theta]$ are all interesting
avenues for future work. Another limitation of our theory is that our proof of
Theorem 2 depends on a conjecture.

We currently only study mixtures with $T$ _independent_ mixture components
without taking into account the cost of producing independent samples of
$\theta$. In reality, this cost depends on the quality of the sampler, warm-up
and burn-in time, and a potentially large number of calls to $\log{\rm
p}(\mathbf{x})$ [Zhang et al., 2021]. Further, $\lambda$ dramatically changes
the shape of $\log\psi(\theta)$, which may affect the efficiency of the
sampler – we mitigated this slightly by scaling the mass parameter of NUTS
with $\lambda$.

We have so far considered $\lambda$ to be constant for a run of our algorithm,
and this can lead to asymptotic bias even when $T$ is large. A simple
adjustment to make our method effective at both small and large $T$ would be
to decay $\lambda$ as $T$ grows, but note that this may require adapting the
sampler parameters on the fly. Our method also requires evaluating ${\rm
KL}({\rm q}||{\rm p})$ many times per sample of $\theta$. This could be made
more efficient by adapting the number of Monte Carlo evaluations (fewer
samples from ${\rm q}$ are sufficient when $\lambda$ is low and components are
narrow), by accounting for stochastic likelihood evaluations [Ma et al.,
2015], or by extending our method to mean-field message-passing [Jaakkola and
Jordan, 1998], where $\nabla_{\theta}{\rm KL}({\rm q}||{\rm p})$ can be
computed in closed form [Hoffman et al., 2013].

###### Acknowledgements.

We thank Emmett Wyman and Roozbeh Farhoudi for helpful discussions early on,
and Konrad Kording for suggestions on writing and presentation. Daniel Lee’s
advice was indispensable for getting our algorithm to run in Stan.

## References

  * Acerbi [2018] Luigi Acerbi.  Variational Bayesian Monte Carlo.  _Advances in Neural Information Processing Systems_ , 2018.
  * Amari [2016] S Amari.  _Information Geometry and Its Applications_.  Applied Mathematical Sciences. Springer Japan, 2016.  ISBN 9784431559788.  URL https://books.google.com/books?id=UkSFCwAAQBAJ.
  * Anaya-Izquierdo and Marriott [2007] Karim Anaya-Izquierdo and Paul Marriott.  Local mixture models of exponential families.  _Bernoulli_ , 13(3):623–640, 2007.  ISSN 1350-7265.  10.3150/07-BEJ6170.  URL http://projecteuclid.org/euclid.bj/1186503479.
  * Angelino et al. [2016] Elaine Angelino, Matthew James Johnson, and Ryan P. Adams.  Patterns of Scalable Bayesian Inference.  _Foundations and Trends in Machine Learning_ , 9(2-3):119–247, 2016.  10.1561/2200000052.
  * Besag et al. [1995] Julian Besag, Peter Green, David Higdon, and Kerrie Mengersen.  Bayesian computation and stochastic systems.  _Statistical Science_ , 10(1):3–44, 1995.  ISSN 08834237.  10.1214/ss/1177010123.
  * Bishop [2006] Christopher M Bishop.  Pattern Recognition and Machine Learning.  _Pattern Recognition_ , page 738, 2006.  ISSN 10179909.  10.1117/1.2819119.  URL http://www.library.wisc.edu/selectedtocs/bg0137.pdf.
  * Blei et al. [2017] David M. Blei, Alp Kucukelbir, and Jon D Mcauliffe.  Variational Inference: A Review for Statisticians.  _arXiv_ , pages 1–41, 2017.
  * Braverman and Bhowmick [2011] Mark Braverman and Abhishek Bhowmick.  Convexity/concavity of mutual information, September 2011.  URL https://www.cs.princeton.edu/courses/archive/fall11/cos597D/L04.pdf.
  * Brunel and Nadal [1998] Nicolas Brunel and Jean Pierre Nadal.  Mutual Information, Fisher Information, and Population Coding.  _Neural Computation_ , 10(7):1731–1757, 1998\.  ISSN 08997667.  10.1162/089976698300017115.
  * Carpenter et al. [2017] Bob Carpenter, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus A. Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.  Stan: A probabilistic programming language.  _Journal of Statistical Software_ , 76(1), 2017.  ISSN 15487660.  10.18637/jss.v076.i01.
  * Chen and Liu [2000] Rong Chen and Jun S Liu.  Mixture kalman filters.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 62(3):493–508, 2000.
  * de Freitas et al. [2001] Nando de Freitas, Pedro Højen-Sørensen, Michael I. Jordan, and Stuart Russel.  Variational MCMC.  _Uncertainty in Artificial Intelligence_ , 2001.
  * Dembo et al. [1991] Amir Dembo, Thomas M. Cover, and Joy A. Thomas.  Information Theoretic Inequalities.  _IEEE Transactions on Information Theory_ , 37(6):1501–1518, 1991.  ISSN 15579654.  10.1109/18.104312.
  * Fauß et al. [2021] Michael Fauß, Alex Dytso, and H. Vincent Poor.  A variational interpretation of the Cramér–Rao bound.  _Signal Processing_ , 182, 2021.  ISSN 01651684.  10.1016/j.sigpro.2020.107917.
  * Gershman et al. [2012] Samuel J. Gershman, Matthew D. Hoffman, and David M. Blei.  Nonparametric Variational Inference.  _Proceedings of the 29th International Conference on Machine Learning_ , pages 235–242, 2012.  ISSN 0899-7667.  10.1162/089976699300016331.  URL https://icml.cc/Conferences/2012/papers/360.pdf.
  * Gill and Levit [1995] Richard D. Gill and Boris Y. Levit.  Applications of the van Trees inequality: A Bayesian Cramér-Rao bound.  _Bernoulli_ , 1(1):59–79, 1995.  URL https://www.jstor.org/stable/3318681.
  * Guo et al. [2016] Fangjian Guo, Xiangyu Wang, Kai Fan, Tamara Broderick, and David B. Dunson.  Boosting Variational Inference.  _arXiv_ , 2016.  URL http://arxiv.org/abs/1611.05559.
  * Harris et al. [2020] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.  10.1038/s41586-020-2649-2.  URL https://doi.org/10.1038/s41586-020-2649-2.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _arXiv preprint arXiv:2006.11239_ , 2020.
  * Hobert and Casella [1996] James P. Hobert and George Casella.  The Effect of Improper Priors on Gibbs Sampling in Hierarchical Linear Mixed Models.  _Journal of the American Statistical Association_ , 91(436):1461–1473, 1996.  ISSN 1537274X.  10.1080/01621459.1996.10476714.
  * Hoffman and Gelman [2014] Matthew D Hoffman and Andrew Gelman.  The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  _Journal of Machine Learning Research_ , 15:1351–1381, 2014\.
  * Hoffman et al. [2013] Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 14:1303–1347, 2013\.  ISSN 1532-4435.  citeulike-article-id:10852147.  URL http://arxiv.org/abs/1206.7051.
  * Hunter [2007] John D. Hunter.  Matplotlib: A 2d graphics environment.  _Computing in Science Engineering_ , 9(3):90–95, 2007.  10.1109/MCSE.2007.55.
  * Jaakkola and Jordan [1998] Tommi S. Jaakkola and Michael I. Jordan.  Improving the Mean Field Approximation via the Use of Mixture Distributions.  In Michael I. Jordan, editor, _Learning in Graphical Models_. Kluwer Academic Publishers, 1998.
  * Kolchinsky and Tracey [2017] Artemy Kolchinsky and Brendan D. Tracey.  Estimating mixture entropy with pairwise distances.  _Entropy_ , 19(7):1–17, 2017.  ISSN 10994300.  10.3390/e19070361.
  * Korattikara et al. [2014] Anoop Korattikara, Yutian Chen, and Max Welling.  Austerity in MCMC Land: Cutting the Metropolis-Hastings Budget.  _International Conference on Machine Learning_ , 32(1):181–189, 2014.  URL http://arxiv.org/abs/1304.5299.
  * Kucukelbir et al. [2017] Alp Kucukelbir, David M. Blei, Andrew Gelman, Rajesh Ranganath, and Dustin Tran.  Automatic Differentiation Variational Inference.  _Journal of Machine Learning Research_ , 18:1–45, 2017\.  ISSN 15337928.
  * Lindsay [1983] Bruce G. Lindsay.  The Geometry of Mixture Likelihoods: A General Theory.  _The Annals of Statistics_ , 11(1):86–94, 1983\.
  * Ma et al. [2015] Yi-An Ma, Tianqi Chen, and Emily B. Fox.  A Complete Recipe for Stochastic Gradient MCMC.  _Advances in Neural Information Processing Systems_ , pages 1–16, 2015.  ISSN 10495258.  URL http://arxiv.org/abs/1506.04696.
  * Magnusson et al. [2021] M. Magnusson, Paul-Christian Bürkner, and Aki Vehtari.  posteriordb: A database of Bayesian posterior inference, 2021.  URL https://github.com/stan-dev/posteriordb.
  * Miller et al. [2017] Andrew C. Miller, Nicholas J. Foti, and Ryan P. Adams.  Variational Boosting: Iteratively Refining Posterior Approximations.  _arXiv_ , 2017.  URL http://arxiv.org/abs/1611.06585.
  * Murphy [2012] Kevin P. Murphy.  _Machine Learning: A Probabilistic Perspective_.  The MIT Press, Cambridge, MA, 2012.
  * Nalisnick and Smyth [2017] Eric Nalisnick and Padhraic Smyth.  Variational Inference with Stein Mixtures.  _NIPS2017 (Workshop)_ , 2017.  ISSN 00368075.  10.1126/science.1070850.  URL https://www.ics.uci.edu/$∼$enalisni/AABI_paper30-Stein_Mixtures.pdf.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Poole et al. [2019] Ben Poole, Sherjil Ozair, Aaron van den Oord, Alexander A. Alemi, and George Tucker.  On variational bounds of mutual information.  _arXiv_ , 2019.  ISSN 23318422.
  * Ranganath et al. [2016] Rajesh Ranganath, Dustin Tran, and David M. Blei.  Hierarchical Variational Models.  _ICML_ , 33:1–9, 2016.
  * Salimans et al. [2015] Tim Salimans, Diederik P. Kingma, and Max Welling.  Markov Chain Monte Carlo and Variational Inference: Bridging the Gap.  _Proceedings of the 32nd International Conference on Machine Learning_ , pages 1218–1226, 2015.  URL http://arxiv.org/abs/1410.6460.
  * Sohl-Dickstein et al. [2015] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _International Conference on Machine Learning_ , pages 2256–2265. PMLR, 2015.
  * Stam [1959] A. J. Stam.  Some inequalities satisfied by the quantities of information of Fisher and Shannon.  _Information and Control_ , 2(2):101–112, 1959\.  ISSN 00199958.  10.1016/S0019-9958(59)90348-1.
  * Stein and Shakarchi [2011] Elias M Stein and Rami Shakarchi.  _Fourier analysis: an introduction_ , volume 1.  Princeton University Press, 2011.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  10.1038/s41592-019-0686-2.
  * Wainwright and Jordan [2008] Martin J. Wainwright and Michael I. Jordan.  Graphical Models, Exponential Families, and Variational Inference.  _Foundations and Trends® in Machine Learning_ , 1(1–2):1–305, 2008.  ISSN 1935-8237.  10.1561/2200000001.
  * Wei and Stocker [2016] Xue-Xin Wei and Alan A. Stocker.  Mutual Information, Fisher Information, and Efficient Coding.  _Neural computation_ , 28(2), 2016.  10.1162/NECO_a_0084.
  * Wolf et al. [2016] Christopher Wolf, Maximilian Karl, and Patrick van der Smagt.  Variational inference with hamiltonian monte carlo.  _arXiv preprint arXiv:1609.08203_ , 2016.
  * Yin and Zhou [2018] Mingzhang Yin and Mingyuan Zhou.  Semi-Implicit Variational Inference.  _International Conference on Machine Learning_ , 35, 2018.
  * Zhang et al. [2019] Cheng Zhang, Judith Butepage, Hedvig Kjellstrom, and Stephan Mandt.  Advances in Variational Inference.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 41(8):2008–2026, 2019.  ISSN 19393539.  10.1109/TPAMI.2018.2889774.
  * Zhang et al. [2021] Lu Zhang, Bob Carpenter, Andrew Gelman, and Aki Vehtari.  Pathfinder: Parallel quasi-Newton variational inference.  _arXiv_ , 2021.  URL http://arxiv.org/abs/2108.03782.
  * Zobay [2014] O. Zobay.  Variational Bayesian inference with Gaussian-mixture approximations.  _Electronic Journal of Statistics_ , 8(1):355–389, 2014.  ISSN 19357524.  10.1214/14-EJS887.

##  A Proofs and Derivations

Throughout, we assume that $\theta$ forms a minimal statistical manifold
[Amari, 2016], so that the degrees of freedom of ${\rm q}$ match the
dimensionality of $\theta$, and whenever ${\rm q}(\mathbf{x};\theta_{i})={\rm
q}(\mathbf{x};\theta_{j})$ for all $\mathbf{x}$, it must be that
$\theta_{i}=\theta_{j}$.

Recall that in the main text, we defined the following objective:

$\mathcal{L}(\psi,\lambda)\equiv\mathcal{I}[\mathbf{x};\theta]-\lambda\mathbb{E}_{\psi(\theta)}\left[{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm p}^{*}(\mathbf{x}))\right]\,,$ ((5)
restated)

where $\lambda\in[1,\infty)$ is a hyper-parameter, and $\psi$ is a probability
density on $\theta$. We also introduced an approximate objective in which
$\mathcal{I}[\mathbf{x};\theta]$ is replaced with

$\mathcal{I}_{\mathcal{F}}[\mathbf{x};\theta]\equiv\mathcal{H}[\theta]-\frac{1}{2}\mathbb{E}_{\psi(\theta)}\left[\log\left|2\pi
e\mathcal{F}(\theta)^{-1}\right|\right]\,.$ ((8) restated)

This approximate objective is

$\mathcal{L}_{\mathcal{F}}(\psi,\lambda)=\mathcal{H}[\theta]+\mathbb{E}_{\psi(\theta)}\left[\frac{1}{2}\log|\mathcal{F}(\theta)|-\lambda{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm p}^{*}(\mathbf{x}))\right]\,,$ ((9)
restated)

and it is maximized for a given $\lambda$ by

$\displaystyle\psi(\theta)$
$\displaystyle=\frac{1}{Z}\exp\left(\frac{1}{2}\log|\mathcal{F}(\theta)|-\lambda{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))\right)$ ((10) restated)
$\displaystyle\text{where}\qquad Z$
$\displaystyle=\int_{\theta}\exp\left(\frac{1}{2}\log|\mathcal{F}(\theta)|-\lambda{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))\right){\rm d}\theta\,.$

###  A.1 Characterizing the Pareto Front

Let us begin with a set of results regarding the shape of the Pareto front
that connects VI to Sampling in Figure 2.

###### Lemma 1

$\mathcal{L}(\psi,\lambda)$ is concave in $\psi$, i.e.
$\mathcal{L}(\omega\psi_{1}+(1-\omega)\psi_{2},\lambda)\geq\omega\mathcal{L}(\psi_{1},\lambda)+(1-\omega)\mathcal{L}(\psi_{2},\lambda)$
for $0\leq\omega\leq 1$. Further, $\mathcal{L}_{\mathcal{F}}(\psi,\lambda)$ is
_strictly_ concave in $\psi$.

#### Proof:

The proof for $\mathcal{L}$ follows from the fact that
$\mathbb{E}_{\psi(\theta)}\left[{\rm KL}({\rm q}(\mathbf{x};\theta)||{\rm
p}^{*}(\mathbf{x}))\right]$ is _linear_ in $\psi$, and
$\mathcal{I}[\mathbf{x};\theta]$ is known to be concave in the marginal
distribution of either variable [Braverman and Bhowmick, 2011]. The proof for
$\mathcal{L}_{\mathcal{F}}$ is similar: the
$\mathbb{E}_{\psi(\theta)}\left[\frac{1}{2}\log|\mathcal{F}(\theta)|\right]$
term is linear in $\psi$, and $\mathcal{H}[\theta]$ is strictly concave in
$\psi$. This can be seen, for instance, by taking the second variational
derivative of $\mathcal{H}[\theta]$ with respect to $\psi$:

$\displaystyle\nabla_{\psi}^{2}\mathcal{H}[\theta]\big{\rvert}_{\theta_{i}\theta_{j}}$
$\displaystyle=\nabla_{\psi}\left(\nabla_{\psi}\mathcal{H}[\theta]\big{\rvert}_{\theta_{i}}\right)\big{\rvert}_{\theta_{j}}$
$\displaystyle=\nabla_{\psi}\left(-\nabla_{\psi}\int_{\theta}\psi(\theta)\log\psi(\theta){\rm
d}\theta\big{\rvert}_{\theta_{i}}\right)\big{\rvert}_{\theta_{j}}$
$\displaystyle=\nabla_{\psi}\left(-1-\log\psi(\theta_{i})\right)\big{\rvert}_{\theta_{j}}$
$\displaystyle=\begin{cases}-\frac{1}{\psi(\theta_{i})}&\text{if
$\theta_{i}=\theta_{j}$}\\\ 0&\text{otherwise}\,.\end{cases}$

Since $\psi(\theta)\geq 0$ everywhere, this implies that the curvature of
$\mathcal{H}[\theta]$ is strictly negative at all values of $\theta$.
$\blacksquare$

###### Lemma 2

Let $\mathcal{I}^{*}(\lambda)$ and $\mathbb{E}[{\rm KL}]^{*}(\lambda)$ denote
the values of Mutual Information and Expected KL achieved by optima of
$\mathcal{L}$ for a given $\lambda$. Then, $\lambda$ defines the slope of the
Pareto front:

$\displaystyle\lambda=\frac{{\rm d}\mathcal{I}^{*}/{\rm d}\lambda}{{\rm
d}\mathbb{E}[{\rm KL}]^{*}/{\rm d}\lambda}\,.$

Or, in the case of $\mathcal{L}_{\mathcal{F}}$, $\lambda$ similarly defines
the slope of

$\displaystyle\lambda=\frac{{\rm d}\mathcal{I}_{\mathcal{F}}^{*}/{\rm
d}\lambda}{{\rm d}\mathbb{E}[{\rm KL}]^{*}/{\rm d}\lambda}\,,$

with $\mathcal{I}_{\mathcal{F}}$ in place of $\mathcal{I}$.

#### Proof:

This follows from viewing $\mathcal{L}$ as the Lagrangian of a constrained
optimization problem, with $\lambda$ as a Lagrange multiplier. The same
argument applies to both $\mathcal{L}$ and $\mathcal{I}$ as to
$\mathcal{L}_{\mathcal{F}}$ and $\mathcal{I}_{\mathcal{F}}$, so we will just
give the proof for one. Consider the constrained optimization problem of
maximizing $\mathcal{I}$ (or $\mathcal{I}_{\mathcal{F}}$) subject to the
constraint that $\mathbb{E}[{\rm KL}({\rm q}||{\rm p})]=C$. The Lagrangian for
this problem is identical to (5), but with $C$ added:

$\displaystyle\mathcal{L}(\psi,\lambda)\equiv\mathcal{I}[\mathbf{x};\theta]-\lambda\left(\mathbb{E}_{\psi(\theta)}\left[{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm p}^{*}(\mathbf{x}))\right]-C\right)$

Optimizing with respect to $\psi$, this is a concave maximization problem with
a linear constraint. A well-known property of such problems is that, at the
solution, the Lagrange multiplier ($\lambda$) is equal to the change in the
objective ($\mathcal{I}^{*}$) per change in the constraint ($C$), or
$\lambda=\frac{{\rm d}\mathcal{I}^{*}}{{\rm d}C}$. Since $C$ is the
constrained value of $\mathbb{E}[{\rm KL}({\rm q}||{\rm p})]$, we also
immediately have $\frac{{\rm d}\mathbb{E}[{\rm KL}]^{*}}{{\rm d}C}=1$. This
implies that

$\displaystyle\lambda=\frac{{\rm d}\mathcal{I}_{\mathcal{F}}^{*}/{\rm
d}C}{{\rm d}\mathbb{E}[{\rm KL}]^{*}/{\rm d}C}\,.$

So far, we have treated $\lambda$ as a function of $C$, but for all values of
$\lambda$ that correspond to a unique $C$, we can invert this relationship and
treat $C$ as a function of $\lambda$. Then, assuming $\frac{{\rm d}C}{{\rm
d}\lambda}\neq 0$ for all $1\leq\lambda<\infty$ that we are interested in, we
have

$\displaystyle\lambda$ $\displaystyle=\frac{{\rm
d}\mathcal{I}_{\mathcal{F}}^{*}/{\rm d}C\times{\rm d}C/{\rm d}\lambda}{{\rm
d}\mathbb{E}[{\rm KL}]^{*}/{\rm d}C\times{\rm d}C/{\rm d}\lambda}=\frac{{\rm
d}\mathcal{I}_{\mathcal{F}}^{*}/{\rm d}\lambda}{{\rm d}\mathbb{E}[{\rm
KL}]^{*}/{\rm d}\lambda}\,.$

Again using the fact that $C=\mathbb{E}[{\rm KL}]^{*}$ by construction, the
condition that $\frac{{\rm d}C}{{\rm d}\lambda}\neq 0$ is equivalent to
$\frac{{\rm d}\mathbb{E}[{\rm KL}]^{*}}{{\rm d}\lambda}\neq 0$. In other
words, as long as changing $\lambda$ has some effect on $\mathbb{E}[{\rm
KL}]^{*}$, the combined effect on $\mathcal{I}^{*}$ and $\mathbb{E}[{\rm
KL}]^{*}$ will be such that $\lambda=\frac{{\rm d}\mathcal{I}^{*}}{{\rm
d}\mathbb{E}[{\rm KL}]^{*}}$. $\blacksquare$

###  A.2 Sampling-like behavior of our method

Recall our definition of sampling:

######  Definition 2 (Sampling)

A stochastic mixture, defined by the component family ${\rm
q}(\mathbf{x};\theta)$ and mixing distribution $\psi(\theta)$, is considered
to be “sampling” if (i) it is unbiased in the limit of infinitely many
components, i.e. ${\rm m}(\mathbf{x})\rightarrow{\rm p}(\mathbf{x})$; and,
(ii) it consists of non-overlapping components. That is, for small values of
of $0<\epsilon\ll 1$, wherever ${\rm q}(\mathbf{x};\theta_{i})>\epsilon$, with
high probability ${\rm q}(\mathbf{x};\theta_{j})<\epsilon$, for all pairs
$\theta_{i},\theta_{j}$ drawn independently from $\psi(\theta)$.

We will assume throughout this section that ${\rm q}$ is a location-scale
family, and in particular Gaussian for Lemma 4, but it may be fruitful for
future work to consider other families of mixture components.

###### Lemma 3

Sampling is an optimum of the original objective, $\mathcal{L}$, when
$\lambda=1$.

#### Proof:

When $\lambda=1$, $\mathcal{L}$ simplifies back to ${\rm KL}({\rm m}||{\rm
p})$. Any unbiased mixture is a minimum of ${\rm KL}({\rm m}||{\rm p})$.
$\blacksquare$

Note, however, that this does not imply sampling is the unique optimum. In
general there may be other unbiased mixing distributions $\psi(\theta)$ such
that ${\rm m}(\mathbf{x})={\rm p}(\mathbf{x})$. For instance, if ${\rm q}$ is
Gaussian and ${\rm p}(\mathbf{x})$ is itself a finite mixture of Gaussians,
then $\psi(\theta)$ could concentrate on exactly those modes in ${\rm p}$. In
any case where there two such unbiased $\psi$s, there are in fact infinitely
many unbiased, since any mixture of them,
$\alpha\psi_{1}(\theta)+(1-\alpha)\psi_{2}(\theta)$, will also be unbiased.
Among all unbiased mixtures, sampling may in some sense be the worst choice –
we conjecture that it has the highest variance of all unbiased mixtures.

###### Lemma 4

When ${\rm q}$ is Gaussian and $\lambda=1$, the optimal $\psi$ that maximizes
the approximate objective $\mathcal{L}_{\mathcal{F}}$ is both unbiased and has
non-overlapping components.

In other words, Lemma 4 states that the solution to the approximate objective
$\mathcal{L}_{\mathcal{F}}$ “looks like” sampling when $\lambda=1$, in the
sense of Definition 1.

#### Proof:

Without loss of generality, let us assume that $\theta$ is already
parameterized in terms of its location and scale,
$[\boldsymbol{\mu},\boldsymbol{\sigma}]$, where $\boldsymbol{\mu}$ determines
the mean of ${\rm q}$ and $\boldsymbol{\sigma}$ determines its covariance.
Then, the Fisher Information Matrix is a block-diagonal
matrix:333https://en.wikipedia.org/wiki/Fisher_information#Multivariate_normal_distribution

$\displaystyle\mathcal{F}(\theta)=\begin{bmatrix}\mathcal{F}(\boldsymbol{\mu})&0\\\
0&\mathcal{F}(\boldsymbol{\sigma})\end{bmatrix}$

where

$\displaystyle\mathcal{F}(\boldsymbol{\mu})$ $\displaystyle=\Lambda$
$\displaystyle\mathcal{F}(\boldsymbol{\sigma})_{ij}$
$\displaystyle=\frac{1}{2}\text{Tr}\left(\Lambda\frac{\partial\Sigma}{\partial\boldsymbol{\sigma}_{i}}\Lambda\frac{\partial\Sigma}{\partial\boldsymbol{\sigma}_{j}}\right)\,.$

$\Lambda$ and $\Sigma$ are the precision matrix and covariance matrix of ${\rm
q}$, respectively. Both $\Lambda$ and $\Sigma$ are functions of the parameters
$\boldsymbol{\sigma}$ but not of $\boldsymbol{\mu}$. To simplify further, let
us assume that the covariance of ${\rm q}$ is diagonal, and that
$\boldsymbol{\sigma}_{i}$ is the log standard deviation of the $i$th dimension
of $\mathbf{x}$:

$\displaystyle\Sigma(\boldsymbol{\sigma})_{ij}=\begin{cases}e^{2\boldsymbol{\sigma}_{i}}&\text{if
$i=j$}\\\ 0&\text{otherwise}\end{cases}$

We emphasize that this simplification is for notational convenience only, and
other parameterizations of $\Sigma(\boldsymbol{\sigma})$ are permissible. With
this assumption, $\mathcal{F}(\boldsymbol{\sigma})$ becomes the identity
matrix, and the log determinant of $\mathcal{F}(\theta)$ becomes simply

$\displaystyle\log|\mathcal{F}(\theta)|=\log|\Lambda|\,.$

So, for Gaussian ${\rm q}$, the expression for $\psi$ becomes

$\displaystyle\log\psi(\theta)=\log\psi(\boldsymbol{\mu},\boldsymbol{\sigma})=\frac{1}{2}\log|\Lambda(\boldsymbol{\sigma})|-\lambda{\rm
KL}({\rm q}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\sigma})||{\rm
p}(\mathbf{x}))\,.$

Next, we will split ${\rm KL}({\rm q}||{\rm p})$ into separate entropy and
cross-entropy terms:

$\displaystyle{\rm KL}({\rm q}||{\rm p})$ $\displaystyle=\mathbb{E}_{{\rm
q}(\mathbf{x};\theta)}\left[\log{\rm
q}(\mathbf{x};\theta)\right]-\mathbb{E}_{{\rm
q}(\mathbf{x};\theta)}\left[\log{\rm p}(\mathbf{x})\right]$
$\displaystyle=-\mathcal{H}[{\rm q}]+\mathcal{CE}[{\rm q}||{\rm p}]\,.$

And note that when ${\rm q}$ is Gaussian, its entropy is given by

$\displaystyle\mathcal{H}[{\rm q}]=\frac{1}{2}\log|2\pi
e\Sigma|=\frac{1}{2}\log|\Sigma|+\text{constants}\,.$

Taking $\lambda=1$ and using the fact that
$\log|\Sigma|=-\log|\Sigma^{-1}|=-\log|\Lambda|$ and combining the above three
equations, the $\mathcal{H}[{\rm q}]$ and
$\log|\mathcal{F}(\boldsymbol{\mu})|$ terms cancel in $\psi$ and we are left –
up to additive constants – with

$\log\psi(\theta)=-\mathcal{CE}[{\rm q}||{\rm p}]=\mathbb{E}_{{\rm
q}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\sigma})}\left[\log{\rm
p}(\mathbf{x})\right]\,.$ (A.1)

To summarize, equation (A.1) says that, using Gaussian components and letting
$\lambda\rightarrow 1$, our method, derived from the
$\mathcal{I}_{\mathcal{F}}$ approximation to $\mathcal{I}$, selects components
simply according to the _cross entropy_ between ${\rm q}(\mathbf{x};\theta)$
and ${\rm p}(\mathbf{x})$.

Note that (A.1) is not a proper distribution over $\theta$. To see this,
consider any sufficiently narrow component such that ${\rm q}$ behaves like a
Dirac delta, or $\mathbb{E}_{{\rm
q}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\sigma})}[\log{\rm
p}(\mathbf{x})]\approx\log{\rm p}(\boldsymbol{\mu})$. Wherever this holds for
some $\boldsymbol{\sigma}$, it will additionally hold for all _narrower_
components at the same $\boldsymbol{\mu}$.444There is an implicit assumption
here that $\log{\rm p}(\mathbf{x})$ is almost everywhere smooth, so that there
is some small enough scale at which ${\rm p}(\mathbf{x})$ appears locally
linear under ${\rm q}$. Therefore, below a particular scale where ${\rm q}$
behaves like a Dirac delta, (A.1) places uniform mass on the infinitely many
${\rm q}$s that are at least as narrow. This effect is visible in the top-
right panel of Figure 2. Also note that $\psi$ is only improper for
$\lambda=1$; for all other $\lambda>1$, a $(\lambda-1)\mathcal{H}[{\rm q}]$
term remains, and $\psi$ cannot place arbitrarily much mass on arbitrarily
narrow components.

Despite its impropriety, we are free to draw samples of $\theta$ from this
improper $\psi$ when $\lambda=1$ [Besag et al., 1995, Hobert and Casella,
1996]. We will then find that with probability approaching $1$ we only ever
see components that “look like” Dirac-deltas. This phenomenon is seen
empirically in all of our experiments where we set $\lambda=1$ and run HMC
dynamics drawing $\theta\sim\psi(\theta)$. Since components will become
arbitrarily narrow, we have the non-overlapping components property required
by our definition of sampling.

Consider decomposing $\psi(\theta)$ into
$\psi(\boldsymbol{\sigma})\psi(\boldsymbol{\mu}|\boldsymbol{\sigma})$. The
previous paragraph establishes that the marginal distribution
$\psi(\boldsymbol{\sigma})$ will allocate effectively all samples to parts of
$\theta$-space where components behave like Dirac deltas. This implies

$\displaystyle\log\psi(\boldsymbol{\mu}|\boldsymbol{\sigma}=\text{narrow})$
$\displaystyle=\mathbb{E}_{{\rm
q}(\mathbf{x};\boldsymbol{\mu},\boldsymbol{\sigma})}\left[\log{\rm
p}(\mathbf{x})\right]$ $\displaystyle=\log{\rm p}(\boldsymbol{\mu})\,.$

Hence, ${\rm m}(\mathbf{x})$ will be a mixture of Dirac-delta-like components,
each of which is chosen in proportion to the true probability of its mean,
${\rm p}(\boldsymbol{\mu})$. This means that ${\rm m}(\mathbf{x})$ will be
unbiased. $\blacksquare$

######  Theorem 3 (Improve on sampling)

If a mixture is sampling as in Definition 1, then $\frac{{\rm d}}{{\rm
d}\lambda}\text{KL bias}=0$ and $\frac{{\rm d}}{{\rm d}\lambda}\text{KL
variance}<0$. Thus, $\frac{{\rm d}}{{\rm d}\lambda}\text{KL error}<0$.

#### Proof:

Our approach will be to calculate the variational derivatives of KL bias and
KL error with respect to $\psi$, then take the inner product (directional
derivative) with the change in $\psi$ per change in $\lambda$.

First, we need the sensitivy of $\psi$ to changes in $\lambda$. Recall that
the closed-form solution for $\psi$ we get from solving
$\mathcal{L}_{\mathcal{F}}$ is

$\displaystyle\log\psi(\theta)=\frac{1}{2}\log|\mathcal{F}(\theta)|-\lambda{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))-\log Z(\lambda)\,.$

The sensitivity of $\log\psi$ to $\lambda$ is

$\displaystyle\frac{{\rm d}}{{\rm d}\lambda}\log\psi(\theta)$
$\displaystyle=-{\rm KL}({\rm q}||{\rm
p})+\frac{1}{Z}\int_{\theta^{\prime}}e^{\frac{1}{2}\log|\mathcal{F}(\theta)|-\lambda{\rm
KL}({\rm q}||{\rm p})}{\rm KL}({\rm q}||{\rm p}){\rm d}\theta^{\prime}$
$\displaystyle=\mathbb{E}_{\psi}[{\rm KL}({\rm q}||{\rm p})]-{\rm KL}({\rm
q}||{\rm p})\,.$

Converting from $\log\psi$ to $\psi$, we get

$\frac{{\rm d}}{{\rm
d}\lambda}\psi(\theta)=\psi(\theta)\left(\mathbb{E}_{\psi}[{\rm KL}({\rm
q}||{\rm p})]-{\rm KL}({\rm q}||{\rm p})\right)$ (A.2)

Recall that we defined $\text{KL bias}={\rm KL}({\rm m}||{\rm p})$ and
$\text{KL variance}=\mathbb{E}[{\rm KL}({\rm m}_{T}||{\rm m})]$. The
variational derivative of Bias with respect to $\psi$, evaluated at
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}$ is

$\displaystyle\nabla_{\psi}{\rm KL}({\rm m}||{\rm p})=$
$\displaystyle\nabla_{\psi}\int_{\mathbf{x}}\left(\mathbb{E}_{\psi}[{\rm
q}(\mathbf{x};\theta)]\right)\log\frac{\left(\mathbb{E}_{\psi}[{\rm
q}(\mathbf{x};\theta)]\right)}{{\rm p}(\mathbf{x})}{\rm d}\mathbf{x}$
$\displaystyle=$ $\displaystyle\int_{\mathbf{x}}\left({\rm
m}(\mathbf{x})\frac{{\rm p}(\mathbf{x})}{{\rm m}(\mathbf{x})}\frac{{\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})}{{\rm
p}(\mathbf{x})}+{\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})\log\frac{{\rm
m}(\mathbf{x})}{{\rm p}(\mathbf{x})}\right){\rm d}\mathbf{x}$ $\displaystyle=$
$\displaystyle 1+\mathbb{E}_{{\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})}\left[\log\frac{{\rm
m}(\mathbf{x})}{{\rm p}(\mathbf{x})}\right]\,.$ (A.3)

To get the sensitivity of Bias to $\lambda$ we will take the inner-product of
(A.2) with (A.3). This is

$\displaystyle\frac{{\rm d}}{{\rm d}\lambda}\text{KL bias}$
$\displaystyle=\left\langle\frac{{\rm d}\text{KL bias}}{{\rm
d}\psi},\frac{{\rm d}\psi}{{\rm d}\lambda}\right\rangle$
$\displaystyle=\int_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}\left(1+\mathbb{E}_{{\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})}\left[\log\frac{{\rm
m}(\mathbf{x})}{{\rm
p}(\mathbf{x})}\right]\right)\psi({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})\left(\mathbb{E}_{\psi}[{\rm
KL}({\rm q}||{\rm p})]-{\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm p})\right){\rm
d}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}$
$\displaystyle=\int_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}(1+0)\psi({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})\left(\mathbb{E}_{\psi}[{\rm
KL}({\rm q}||{\rm p})]-{\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm p})\right){\rm
d}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}$
(unbiased) $\displaystyle=\mathbb{E}_{\psi}[{\rm KL}({\rm q}||{\rm
p})]-\mathbb{E}_{\psi}[{\rm KL}({\rm q}||{\rm p})]$ $\displaystyle=0\,.$

So, we can conclude that in the sampling limit, small changes in $\lambda$
have no effect on Bias. Geometrically, this tells us the Pareto Front is
tangent to the y=x line in that limit, as illustrated in Figure 2.

Next we will consider the variational derivative of the Variance component of
KL error with respect to $\psi$, where

KL variance $\displaystyle\equiv\mathbb{E}_{1..T}[{\rm KL}({\rm m}_{T}||{\rm
m})]$
$\displaystyle=\mathbb{E}_{1..T}\left[\int_{\mathbf{x}}\left(\frac{1}{T}\sum_{t=1}^{T}{\rm
q}(\mathbf{x};\theta_{t})\right)\log\frac{\left(\frac{1}{T}\sum_{j=1}^{T}{\rm
q}(\mathbf{x};\theta_{j})\right)}{{\rm m}(\mathbf{x})}{\rm
d}\mathbf{x}\right]$

using the shorthand $\mathbb{E}_{1..T}[\ldots]$ to denote an expectation over
independent draws of $\\{{\theta_{t}}\\}\sim\psi(\theta)$. We will apply the
assumption of non-overlapping components to simplify ${\rm KL}({\rm
m}_{T}||{\rm m})$. Let $\int_{\mathbf{x}\in{\rm q}_{t}}\ldots{\rm
d}\mathbf{x}$ denote an integral over just the region of $\mathbf{x}$-space
where ${\rm q}(\mathbf{x};\theta_{t})>\epsilon$ for some small $\epsilon$. By
assumption, these regions are disjoint for all pairs of $\theta$s, with high
probability. Splitting the integral into $T$ separate regions and rearranging
terms inside the $\log$, we have

KL variance
$\displaystyle\approx\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}_{1..T}\left[\int_{\mathbf{x}\in{\rm
q}_{t}}{\rm q}(\mathbf{x};\theta_{t})\log\left(\frac{{\rm
q}(\mathbf{x};\theta_{t})}{{\rm
m}(\mathbf{x})}\left(\frac{1}{T}+\frac{1}{T}\sum_{j\neq t}\frac{{\rm
q}(\mathbf{x};\theta_{j})}{{\rm q}(\mathbf{x};\theta_{t})}\right)\right){\rm
d}\mathbf{x}\right]+\mathcal{O}(\epsilon)$
$\displaystyle=\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}_{1..T}\left[\int_{\mathbf{x}\in{\rm
q}_{t}}{\rm q}(\mathbf{x};\theta_{t})\left(\log\frac{{\rm
q}(\mathbf{x};\theta_{t})}{{\rm
m}(\mathbf{x})}+\log\left(\frac{1}{T}+\frac{1}{T}\sum_{j\neq t}\frac{{\rm
q}(\mathbf{x};\theta_{j})}{{\rm q}(\mathbf{x};\theta_{t})}\right)\right){\rm
d}\mathbf{x}\right]+\mathcal{O}(\epsilon)$
$\displaystyle=\mathbb{E}_{\psi}[{\rm KL}({\rm q}||{\rm
m})]+\frac{1}{T}\sum_{t=1}^{T}\mathbb{E}_{1..T}\left[\int_{\mathbf{x}\in{\rm
q}_{t}}{\rm
q}(\mathbf{x};\theta_{t})\log\left(\frac{1}{T}+\frac{1}{T}\sum_{j\neq
t}\frac{{\rm q}(\mathbf{x};\theta_{j})}{{\rm
q}(\mathbf{x};\theta_{t})}\right){\rm
d}\mathbf{x}\right]+\mathcal{O}(\epsilon)$

From here, we can get an upper-bound on Variance by noting that $\frac{{\rm
q}(\mathbf{x};\theta_{j})}{{\rm q}(\mathbf{x};\theta_{t})}\leq 1$ by the non-
overlapping assumption. Since there are $T-1$ of these in the sum,
$\log\left(\frac{1}{T}+\frac{1}{T}\sum_{j\neq t}\frac{{\rm
q}(\mathbf{x};\theta_{j})}{{\rm
q}(\mathbf{x};\theta_{t})}\right)\leq\left(\frac{1}{T}+\frac{T-1}{T}\right)=0$,
and so

$\displaystyle\text{KL variance}\leq\mathbb{E}_{\psi}[{\rm KL}({\rm q}||{\rm
m})]+\mathcal{O}(\epsilon)\,.$

Note that this bound used the assumption of non-overlapping components, and is
therefore only applicable for small $\lambda$ and moderate values of $T$.
Since we are interested in showing that $\frac{{\rm d}}{{\rm
d}\lambda}\text{KL variance}<0$ in the sampling limit, showing that the _upper
bound_ on variance decreases with $\lambda$ will suffice. Using this upper-
bound, we get the following variational derivative of Variance with respect to
$\psi$ at each value of
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}$:

$\displaystyle\nabla_{\psi}\text{KL
variance}\big{\rvert}_{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}}$
$\displaystyle\approx\nabla_{\psi}\left.\int_{\theta}\psi(\theta)\int_{\mathbf{x}}{\rm
q}(\mathbf{x};\theta)\log\frac{{\rm q}(\mathbf{x};\theta)}{{\rm
m}(\mathbf{x})}{\rm d}\mathbf{x}{\rm
d}\theta\right\rvert_{{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}}$
$\displaystyle=-\int_{\theta}\psi(\theta)\int_{\mathbf{x}}{\rm
q}(\mathbf{x};\theta)\frac{{\rm m}(\mathbf{x})}{{\rm
q}(\mathbf{x};\theta)}\frac{{\rm q}(\mathbf{x};\theta)}{{\rm
m}(\mathbf{x})^{2}}{\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}){\rm
d}\mathbf{x}{\rm d}\theta+\int_{\mathbf{x}}{\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})\log\frac{{\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})}{{\rm
m}(\mathbf{x})}{\rm d}\mathbf{x}$ $\displaystyle=-1+{\rm KL}({\rm
q}(\mathbf{x};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})||{\rm
m}(\mathbf{x}))\,.$

Taking the inner product with $\frac{{\rm d}}{{\rm d}\lambda}\psi$, and
applying the unbiased assumption,

$\displaystyle\frac{{\rm d}}{{\rm d}\lambda}\text{KL variance}$
$\displaystyle=\left\langle\frac{{\rm d}\text{KL variance}}{{\rm
d}\psi},\frac{{\rm d}\psi}{{\rm d}\lambda}\right\rangle$
$\displaystyle=\int_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}\left(-1+{\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm
m})\right)\psi({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})\left(\mathbb{E}_{\psi}[{\rm
KL}({\rm q}||{\rm p})]-{\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm p})\right){\rm
d}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}$
$\displaystyle=\int_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}\left(-1+{\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm
p})\right)\psi({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})\left(\mathbb{E}_{\psi}[{\rm
KL}({\rm q}||{\rm p})]-{\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm p})\right){\rm
d}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta}$
(unbiased)
$\displaystyle=-\mathbb{E}_{\psi({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\theta})}\left[\left({\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm p})-\mathbb{E}_{\psi}[{\rm KL}({\rm q}||{\rm p})]\right){\rm
KL}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}{\rm
q}}||{\rm p})\right]$ $\displaystyle=-\text{var}\left({\rm KL}({\rm q}||{\rm
p})\right)\,.$

In other words, this says that the change in the (upper bound on) “Variance,”
defined as $\mathbb{E}_{1..T}{\rm KL}({\rm m}_{T}||{\rm m})]$, is _negative_ ,
with magnitude given by the variance of the values taken by ${\rm KL}({\rm
q}||{\rm p})$ across all $\theta$.

To summarize, we have shown that, in the sampling limit, where $\lambda=1$, we
have $\frac{{\rm d}}{{\rm d}\lambda}\text{KL bias}=0$ and $\frac{{\rm d}}{{\rm
d}\lambda}\text{KL variance}\leq 0$, which proves the lemma. $\blacksquare$

###  A.3 VI-like behavior of our method

######  Definition 3 (VI limit)

We model the large $\lambda$ limit of our method using a Laplace approximation
around the optimal $\theta^{*}=\operatorname*{arg\,min}_{\theta}{\rm KL}({\rm
q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))$:

$\begin{split}\psi(\theta)&\approx\mathcal{N}(\theta;\theta^{*},{\Sigma^{*}})\\\
\text{where}\qquad{\Sigma^{*}}^{-1}&=\lambda\nabla^{2}_{\theta}{\rm KL}({\rm
q}(\mathbf{x};\theta)||{\rm
p}(\mathbf{x}))\big{\rvert}_{\theta^{*}}\,.\end{split}$ (A.4)

In other words, we approximate $\psi$ by a normal distribution whose mean is
$\theta^{*}$ and whose precision is set by the curvature of ${\rm KL}({\rm
q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))$ and scales with $\lambda$. We
will assume, for the purposes of proofs related to the VI limit, that there is
a single optimal $\theta^{*}$.

######  Theorem 4 (Improve on VI)

Assume that ${\rm p}(\mathbf{x})$ is heaver-tailed than ${\rm
q}(\mathbf{x};\theta^{*})$. Then, there exists some $T_{0}>1$ such that for
all $T\geq T_{0}$, $\frac{{\rm d}}{{\rm d}\lambda}\text{KL error}>0$, when
$\lambda$ is sufficiently large.

#### Proof:

As $\lambda$ grows, the Laplace approximation in (A.4) becomes increasingly
narrow. This allows us to approximate expectations under $\psi$ using a second
order Taylor approximation to the integrand. The general rule for multivariate
Gaussians is

$\displaystyle\mathbb{E}_{\mathcal{N}(\mathbf{y};\mu,\Sigma)}[f(\mathbf{y})]\approx
f(\mu)+\frac{1}{2}\text{Tr}\left(\Sigma\;\nabla^{2}_{\mathbf{y}}f\right)\big{\rvert}_{\mu}$

Recall that we defined KL error as $\mathbb{E}_{1..T}[{\rm KL}({\rm
m}_{T}(\mathbf{x})||{\rm p}(\mathbf{x}))]$. Approximating each
$\psi(\theta_{t})$ as a multivariate Gaussian, their product is also a
multivariate Gaussian whose collective covariance is block-diagonal555This
assumes the $T$ components are statistically independent draws from
$\psi(\theta)$. The approach outlined here could be generalized to include
correlation between $\theta$s in the off-block-diagonals to model variance of
an autocorrelated chain of $\theta$ values. containing $T$ copies of
${\Sigma^{*}}$ from (A.4), and whose collective mean is $\theta^{*}$ for each
component. At this mean value where all $T$ components’ parameters are equal
to $\theta^{*}$, ${\rm m}_{T}(\mathbf{x})$ becomes ${\rm
q}(\mathbf{x};\theta^{*})$. Hence, applying the Taylor series approximation to
KL error, the $f(\mu)$ term is just ${\rm KL}({\rm
q}(\mathbf{x};\theta^{*})||{\rm p}(\mathbf{x}))$. The second term is

$\displaystyle\frac{1}{2}\text{Tr}\left(\begin{bmatrix}{\Sigma^{*}}&&&0\\\
&{\Sigma^{*}}&&&\\\ &&\ddots&\\\
0&&&{\Sigma^{*}}\end{bmatrix}\nabla^{2}_{\theta_{1},\ldots,\theta_{T}}{\rm
KL}({\rm m}_{T}||{\rm p})\right)\,.$

First, note that the zeros in the off-block-diagonal terms on the left mean
that we can ignore interactions between $\theta$s across different mixture
components in the Hessian term on the right. Second, there is $T-$fold
symmetry between all components. So, this simplifies to

$\displaystyle\frac{T}{2}\text{Tr}\left({\Sigma^{*}}\;\nabla^{2}_{\theta_{1}}{\rm
KL}({\rm m}_{T}||{\rm
p})\right)=\frac{T}{2\lambda}\text{Tr}\left((\nabla^{2}_{\theta}{\rm KL}({\rm
q}||{\rm p}))^{-1}\;\nabla^{2}_{\theta_{1}}{\rm KL}({\rm m}_{T}||{\rm
p})\right)\,.$

Next, since this Hessian is being evaluated around the point $\theta^{*}$, all
of $\theta_{2},\ldots,\theta_{T}$ are equal to $\theta^{*}$, and we can write
the mixture as a function only of the component parameters we are varying in
the Hessian. Call this mixture with $T-1$ components set to the variational
solution ${\rm m}_{T}^{*}$, defined as

$\displaystyle{\rm m}_{T}^{*}(\mathbf{x};\theta)=\frac{T-1}{T}{\rm
q}(\mathbf{x};\theta^{*})+\frac{1}{T}{\rm q}(\mathbf{x};\theta)\,.$

We will now calculate each of these Hessians. Note: in what follows we will
use $\theta_{i}$ and $\theta_{j}$ to indicate the $i$ and $j$th indices of the
vector $\theta$, whereas we had used $\theta_{t}$ to indicate one of $T$
vectors. For ${\Sigma^{*}}$, we need the second derivative (Hessian) of ${\rm
KL}({\rm q}||{\rm p})$:

$\displaystyle\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{i}}{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))$
$\displaystyle=\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{i}}\int_{\mathbf{x}}{\rm
q}(\mathbf{x};\theta)\log\frac{{\rm q}(\mathbf{x};\theta)}{{\rm
p}(\mathbf{x})}{\rm d}\mathbf{x}$
$\displaystyle=\frac{\partial}{\partial\theta_{j}}\int_{\mathbf{x}}\left[\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)\left(1+\log\frac{{\rm q}(\mathbf{x};\theta)}{{\rm
p}(\mathbf{x})}\right)\right]{\rm d}\mathbf{x}$
$\displaystyle=\int_{\mathbf{x}}\left[\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)\left(\frac{\frac{\partial}{\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)}{{\rm
q}(\mathbf{x};\theta^{*})}\right)+\left(\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)\right)\left(1+\log\frac{{\rm
q}(\mathbf{x};\theta^{*})}{{\rm p}(\mathbf{x})}\right)\right]{\rm
d}\mathbf{x}$ $\displaystyle(*)$
$\displaystyle=\int_{\mathbf{x}}\frac{\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)\left(\frac{\partial}{\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)\right)}{{\rm q}(\mathbf{x};\theta^{*})}{\rm
d}\mathbf{x}+\int_{\mathbf{x}}\left(\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)\right)\log\frac{{\rm q}(\mathbf{x};\theta^{*})}{{\rm
p}(\mathbf{x})}{\rm d}\mathbf{x}$
$\displaystyle=\mathcal{F}(\theta^{*})+M(\theta^{*})\,.$ (A.5)

In line $(*)$ we used the fact that $\int_{\mathbf{x}}\nabla^{2}_{\theta}{\rm
q}(\mathbf{x};\theta){\rm
d}\mathbf{x}=\nabla^{2}_{\theta}\int_{\mathbf{x}}{\rm
q}(\mathbf{x};\theta){\rm d}\mathbf{x}=\nabla^{2}_{\theta}1=0$. $\mathcal{F}$
is the Fisher Information Matrix, and we have defined
$M(\theta)=\int_{\mathbf{x}}\left(\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)\right)\log\frac{{\rm q}(\mathbf{x};\theta)}{{\rm
p}(\mathbf{x})}{\rm d}\mathbf{x}$.

Following a similar derivation, the Hessian of ${\rm KL}({\rm
m}_{T}^{*}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))$ is

$\displaystyle\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{i}}$
$\displaystyle{\rm KL}({\rm m}_{T}^{*}(\mathbf{x};\theta)||{\rm
p}(\mathbf{x}))=\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{i}}\int_{\mathbf{x}}\left(\frac{T-1}{T}{\rm
q}(\mathbf{x};\theta^{*})+\frac{1}{T}{\rm
q}(\mathbf{x};\theta)\right)\log\frac{\left(\frac{T-1}{T}{\rm
q}(\mathbf{x};\theta^{*})+\frac{1}{T}{\rm q}(\mathbf{x};\theta)\right)}{{\rm
p}(\mathbf{x})}{\rm d}\mathbf{x}$
$\displaystyle=\frac{\partial}{\partial\theta_{j}}\int_{\mathbf{x}}\left[\frac{1}{T}\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)+\frac{1}{T}\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)\log\frac{\left(\frac{T-1}{T}{\rm
q}(\mathbf{x};\theta^{*})+\frac{1}{T}{\rm q}(\mathbf{x};\theta)\right)}{{\rm
p}(\mathbf{x})}\right]{\rm d}\mathbf{x}$
$\displaystyle=\frac{1}{T}\frac{\partial}{\partial\theta_{j}}\int_{\mathbf{x}}\left[\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)\left(1+\log\frac{\left(\frac{T-1}{T}{\rm
q}(\mathbf{x};\theta^{*})+\frac{1}{T}{\rm q}(\mathbf{x};\theta)\right)}{{\rm
p}(\mathbf{x})}\right)\right]{\rm d}\mathbf{x}$
$\displaystyle=\frac{1}{T}\int_{\mathbf{x}}\left[\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)\left(\frac{\frac{1}{T}\frac{\partial}{\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)}{{\rm
m}_{T}^{*}(\mathbf{x};\theta)}\right)+\left(\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)\right)\left(1+\log\frac{{\rm
m}_{T}^{*}(\mathbf{x};\theta)}{{\rm p}(\mathbf{x})}\right)\right]{\rm
d}\mathbf{x}$ $\displaystyle(**)$
$\displaystyle=\frac{1}{T^{2}}\int_{\mathbf{x}}\frac{\left(\frac{\partial}{\partial\theta_{i}}{\rm
q}(\mathbf{x};\theta)\right)\left(\frac{\partial}{\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)\right)}{{\rm q}(\mathbf{x};\theta^{*})}{\rm
d}\mathbf{x}+\frac{1}{T}\int_{\mathbf{x}}\left(\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}{\rm
q}(\mathbf{x};\theta)\right)\log\frac{{\rm q}(\mathbf{x};\theta^{*})}{{\rm
p}(\mathbf{x})}{\rm d}\mathbf{x}$
$\displaystyle=\frac{1}{T^{2}}\mathcal{F}(\theta^{*})+\frac{1}{T}M(\theta^{*})$
$\displaystyle=\frac{1}{T}\frac{\partial^{2}}{\partial\theta_{j}\partial\theta_{i}}{\rm
KL}({\rm q}(\mathbf{x};\theta)||{\rm
p}(\mathbf{x}))+\mathcal{F}(\theta)\left(\frac{1-T}{T^{2}}\right)$ (A.6)

Here, in $(**)$, we additionally used the fact that ${\rm
m}_{T}^{*}(\mathbf{x};\theta^{*})={\rm q}(\mathbf{x};\theta^{*})$. We then
wrote the final line in terms of the Hessian of ${\rm KL}({\rm q}||{\rm p})$
in (A.5).

To summarize, near the variational limit we have that the KL error is
approximately

$\displaystyle{\rm KL}({\rm q}(\mathbf{x};\theta^{*})||{\rm
p}(\mathbf{x}))+\frac{T}{2\lambda}\text{Tr}((\underbrace{\nabla^{2}_{\theta}{\rm
KL}({\rm q}||{\rm
p})}_{(\ref{eqn:hess_kl_q_p})})^{-1}\;(\underbrace{\nabla^{2}_{\theta}{\rm
KL}({\rm m}_{T}^{*}||{\rm p})}_{(\ref{eqn:hess_kl_mt_p})}))\,.$

Plugging in (A.5) and (A.6), this is

KL error $\displaystyle\approx{\rm KL}({\rm q}(\mathbf{x};\theta^{*})||{\rm
p}(\mathbf{x}))+\frac{1}{2\lambda}\text{Tr}\left(\mathbf{I}+\frac{1-T}{T}(\mathcal{F}+M)^{-1}\mathcal{F}\right)$
$\displaystyle={\rm KL}({\rm q}(\mathbf{x};\theta^{*})||{\rm
p}(\mathbf{x}))+\frac{d}{2\lambda}-\frac{1}{2\lambda}\text{Tr}\left(\frac{T-1}{T}(\mathcal{F}+M)^{-1}\mathcal{F}\right)$

where $\mathbf{I}$ is the identity matrix. Consider the case where $T=1$: the
KL error simplifies to ${\rm KL}({\rm q}(\mathbf{x};\theta^{*})||{\rm
p}(\mathbf{x}))+\frac{d}{2\lambda}$ where $d$ is the dimensionality of
$\theta$. Therefore when $T=1$, KL error is only reduced by further increasing
$\lambda$. This is an intuitive result: we cannot reduce bias compared to VI
when using a single component, and any stochasticity only adds variance.

Now consider the case where $T\geq 2$. We are interested in cases where KL
error _increases_ with $\lambda$ near the VI limit. This is equivalent to
asking when the following inequality holds:

$\displaystyle\text{Tr}\left((\mathcal{F}+M)^{-1}\mathcal{F}\right)>\frac{T}{T-1}\overbrace{\text{Tr}(\mathbf{I})}^{d}\,,$

Recall from (A.5) that $\mathcal{F}+M$ is the Hessian of ${\rm KL}({\rm
q}(\mathbf{x};\theta)||{\rm p}(\mathbf{x}))$, and note that the Fisher
Information Matrix is equivalent to the Hessian with respect to $\theta$ of
${\rm KL}({\rm q}(\mathbf{x};\theta)||{\rm q}(\mathbf{x};\theta^{*}))$, so we
can rewrite this inequality as

$\displaystyle\text{Tr}\left((\nabla^{2}_{\theta}{\rm KL}({\rm q}||{\rm
p}))^{-1}\nabla^{2}_{\theta}{\rm KL}({\rm q}||{\rm
q}^{*})\right)>\frac{T}{T-1}\text{Tr}(\mathbf{I})\,.$

Assuming $\theta^{*}$ is a local minimum of ${\rm KL}({\rm q}||{\rm p})$
(which follows from the assumption that $\theta^{*}$ is the unique minimum),
both of these are positive definite matrices encoding how sharply curved the
${\rm KL}({\rm q}||{\rm p})$ or ${\rm KL}({\rm q}||{\rm q}^{*})$ objectives
are.

If ${\rm p}$ is in the same family as ${\rm q}$, then ${\rm q}*={\rm p}$ and
this inequality becomes
$\text{Tr}\left(\mathbf{I}\right)>\frac{T}{T-1}\text{Tr}(\mathbf{I})$, which
is false for all finite $T$ and approaches equality as $T\rightarrow\infty$.
This again captures the intuitive idea that we cannot improve on VI by
reducing $\lambda$ when the single-component ${\rm q}$ is already unbiased.

Conversely, we can view the ratio

$\displaystyle\frac{\text{Tr}\left((\nabla^{2}_{\theta}{\rm KL}({\rm q}||{\rm
p}))^{-1}\nabla^{2}_{\theta}{\rm KL}({\rm q}||{\rm
q}^{*})\right)}{\text{Tr}(\mathbf{I})}$ (A.7)

as an indication of _how poorly matched_ ${\rm q}(\mathbf{x};\theta^{*})$ is
to ${\rm p}(\mathbf{x})$, locally around the single-component VI solution. We
conjecture that this ratio is always greater than $1$ whenever ${\rm
p}(\mathbf{x})$ is heavier-tailed than ${\rm q}(\mathbf{x};\theta^{*})$. Since
$\frac{T}{T-1}$ approaches $1$ from above in the limit of large $T$, this
implies that there will be some finite $T_{0}$ where $\frac{T_{0}}{T_{0}-1}$
is less than the ratio in (A.7), and that such a $T_{0}$ will be reached
sooner the worse ${\rm q}(\mathbf{x};\theta^{*})$ locally approximates ${\rm
p}$. $\blacksquare$

##  B Additional Experiments

Figure B.1: Further examples. To verify that the results in the main
manuscript generalize other datasets, we ran our algorithm at multiple values
of $\lambda$ on three problems from the posteriordb dataset [Magnusson et al.,
2021]. Rather than show convergence for one random $f(\mathbf{x})$, as in the
main text, we generated 200 random functions over the unconstrained parameter
space of each model, using the Fourier synthesis method described in section C
below, with $\alpha=-1$. We then calculated the expectation $\mathbb{E}_{{\rm
m}_{T}(\mathbf{x})}[f(\mathbf{x})]$ after subsampling mixtures of $T=100$
components, and plotted the distribution of squared error relative to ground
truth expectations (based on very long runs of Stan’s default NUTS
implementation). Error bars reflect the combined effect of two sources of
variability: one from the random choice of $f(\mathbf{x})$ and one from the
random subsampling of mixtures. Blue line is the median squared error across
different $f$s and different mixtures, and shading shows $[25\%,75\%]$
quantiles. The blue dot to the right represents the expectation calculated by
a mean-field VI approximation (a diagonal Gaussian) using the automatic
differentiation variational inference (ADVI) package built in to Stan
[Kucukelbir et al., 2017], run with its default parameters.

##  C Numerical Details

We implemented (10) in Stan [Carpenter et al., 2017]. For ${\rm q}$, we used
the family of multivariate Gaussians with diagonal covariance, parameterized
as $\theta=[\mu_{1},\ldots,\mu_{n},\log\sigma_{1},\ldots,\log\sigma_{n}]$
where $n$ is the number of unconstrained parameters (i.e the dimensionality of
$\mathbf{x}$). In this parameterization, $\frac{1}{2}\log\mathcal{F}(\theta)$
is simply $-\sum_{i=1}^{n}\log\sigma_{i}$. We sampled $\theta$ from
$\psi(\theta)$ using Stan’s default implementation of the No U-Turn Sampler
(NUTS) with automatic step-size adaptation [Hoffman and Gelman, 2014], and we
set the mass equal to $\lambda$ times the identity matrix. NUTS requires both
${\rm KL}({\rm q}||{\rm p})$ and its gradient, which we computed using Monte
Carlo samples from ${\rm q}$ and the reparameterization trick. The
reparameterized samples were frozen for each trajectory of NUTS and resampled
between trajectories.

All code to generate the figures in this paper is available publicly online;
the repository URL will be shared after the double-blind review process is
complete. Python libraries used include NumPy, SciPy, PyTorch, and Matplotlib
[Harris et al., 2020, Virtanen et al., 2020, Paszke et al., 2019, Hunter,
2007].

###  C.1 Figure details

We used two toy distributions in our results:

  * •

The “banana” distribution over $\mathbb{R}^{2}$, defined as

$\displaystyle\log{\rm p}(x,y)=-(y-(x/2)^{2})^{2}-(x/2)^{2}\,.$

  * •

The “Laplace mixture” distribution over $\mathbb{R}^{1}$, defined as

$\displaystyle{\rm p}(x)\propto
0.4e^{\frac{|x+1.5|}{0.75}}+0.6e^{\frac{|x-1.5|}{0.75}}\,.$

We also tested our method on three reference problems taken from posteriordb
[Magnusson et al., 2021], a database of reference problems for testing and
validating inference methods. These were arK, eigh schools centered, and
garch11. Results for these additional problems are shown in Figure B.1.

In our experiments, all functions integrated are sums of sinusoids,

$f(\mathbf{x})=\sum_{\omega=1}^{N}a\sin(\omega\mathbf{t}^{T}\mathbf{x}+\phi_{\omega})$

where $\mathbf{t}$ is a random unit vector. This is a convenient target
distribution as the integral of a sinusoid under a Gaussian is known
analytically:

$\int_{\mathbf{x}}\sin(\omega\mathbf{t}^{\top}\mathbf{x}+\phi_{\omega})\mathcal{N}(\mu,\Sigma)=\sin(\omega\mathbf{t}^{\top}\mu+\phi_{\omega})\exp\left(-\frac{\omega^{2}}{2}\mathbf{t}^{T}\Sigma\mathbf{t}\right)$

The capability for exact integration of $\int_{\mathbf{x}}{\rm
m}_{T}(\mathbf{x})f(\mathbf{x}){\rm d}\mathbf{x}$ ensures that no additional
variance is introduced in plots; all variance is due to the selection of the
components ${\rm q}$. In general this integral can be computed with MC methods
or, in low enough dimensions, Gaussian quadrature.

In our experiments (Figures 3 and 4) we used $N=100$ sinusoidal components in
$f(\mathbf{x})$, and calculated bias using $T=5,000$ components thinned from 4
MCMC chains of length $50,000$. To calculate variance, we subsampled $T=10$
components from these chains, and computed variance over these random
instantiations of ${\rm m}_{10}(\mathbf{x})$. The NUTS samples over
$\mathbf{x}$ treated as ground truth derive from 4 chains of length 1,000,000.